15 — Géométrie Computationnelle : Cours complet

Niveau : Université — Durée : 4h de cours + 4h de TP La géométrie computationnelle est au cœur des SIG, du graphisme, de la robotique et des jeux vidéo.

---

Partie I — Points et vecteurs

1.1 Représentation

Point  = (x, y)
Vecteur = (dx, dy) = B − A

Un point et un vecteur sont stockés de la même façon ; le sens diffère.

def add(a, b): return (a[0] + b[0], a[1] + b[1])
def sub(a, b): return (a[0] - b[0], a[1] - b[1])
def scale(v, s): return (v[0] * s, v[1] * s)

1.2 Distance et norme

norm(v) = sqrt(vx² + vy²)
dist(A, B) = norm(B − A)

---

Partie II — Produits scalaire et vectoriel

2.1 Produit scalaire (dot product)

a · b = ax·bx + ay·by = |a|·|b|·cos(θ)

• Signe du produit scalaire → angle :

  • 0 : angle aigu

  • = 0 : perpendiculaires
  • < 0 : angle obtus

def dot(a, b):
    return a[0] * b[0] + a[1] * b[1]

2.2 Produit vectoriel (cross product) — 2D

cross(a, b) = ax·by − ay·bx = |a|·|b|·sin(θ)

• Signe du cross → sens de rotation :

  • 0 : b est à gauche de a (rotation antihoraire)

  • < 0 : b est à droite de a (rotation horaire)
  • = 0 : colinéaires

def cross(a, b):
    return a[0] * b[1] - a[1] * b[0]

Le produit vectoriel 2D est la clé de voûte de toute la géométrie computationnelle : orientation, intersection, aire, convex hull.

---

Partie III — Orientation de trois points

3.1 Définition

L'orientation de (a, b, c) est le signe de cross(b − a, c − a) :

def orient(a, b, c):
    return cross(sub(b, a), sub(c, a))

• > 0 : tourne à gauche → CCW (counterclockwise)
• < 0 : tourne à droite → CW (clockwise)
• = 0 : les 3 points sont collinéaires

Diagramme en cours de génération...

3.2 Précision numérique

⚠️ Les calculs en virgule flottante sont imprécis. Deux stratégies :

1. Entiers quand possible (coordonnées entières → cross exact en entiers).
2. Epsilon : |orient| < EPS ⇒ collinéaires, avec EPS ~ 1e-9.

---

Partie IV — Segments

4.1 Point sur segment

c est sur le segment [a, b] ssi :

1. orient(a, b, c) == 0 (collinéaire) ;
2. c dans la boîte englobante : min(ax,bx) ≤ cx ≤ max(ax,bx) et idem pour y.

def on_segment(a, b, c):
    if orient(a, b, c) != 0:
        return False
    return (min(a[0], b[0]) <= c[0] <= max(a[0], b[0]) and
            min(a[1], b[1]) <= c[1] <= max(a[1], b[1]))

4.2 Intersection de segments

Deux segments s'intersectent si :

• Cas général : les orientations sont opposées : orient(a,b,c)·orient(a,b,d) < 0 ET orient(c,d,a)·orient(c,d,b) < 0.
• Cas dégénérés (collinéarité) : un point est sur l'autre segment.

def segments_intersect(a, b, c, d):
    o1 = orient(a, b, c)
    o2 = orient(a, b, d)
    o3 = orient(c, d, a)
    o4 = orient(c, d, b)
    if ((o1 > 0) != (o2 > 0)) and ((o3 > 0) != (o4 > 0)):
        return True                       # cas général strict
    if o1 == 0 and on_segment(a, b, c): return True   # dégénérés
    if o2 == 0 and on_segment(a, b, d): return True
    if o3 == 0 and on_segment(c, d, a): return True
    if o4 == 0 and on_segment(c, d, b): return True
    return False

---

Partie V — Polygones

5.1 Aire d'un polygone (formule du lacet / shoelace)

def polygon_area(points):
    s = 0
    n = len(points)
    for i in range(n):
        x1, y1 = points[i]
        x2, y2 = points[(i + 1) % n]
        s += x1 * y2 - x2 * y1
    return abs(s) / 2

Signe : si les points sont en ordre CCW, s > 0 (aire orientée positive) ; CW → négatif.

5.2 Point in polygon (ray casting)

Lancer un rayon vers +∞ ; compter les intersections avec les arêtes. Impair = à l'intérieur, pair = à l'extérieur.

def point_in_polygon(pt, poly):
    x, y = pt
    inside = False
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]
        x2, y2 = poly[(i + 1) % n]
        if ((y1 > y) != (y2 > y)) and (x < (x2 - x1) * (y - y1) / (y2 - y1) + x1):
            inside = not inside
    return inside

5.3 Pour les polygones convexes

Alternative O(log n) : bissecter par rapport au premier sommet et vérifier l'orientation (binary search).

---

Partie VI — Convex Hull

6.1 Définition

Le convex hull d'un ensemble de points est le plus petit polygone convexe les contenant tous.

Diagramme en cours de génération...

6.2 Graham scan (O(n log n))

1. Trier par angle autour du point le plus bas.
2. Pour chaque point, tant que le dernier virage n'est pas à gauche (CCW), pop du stack.
3. Pousser le point.

def convex_hull_graham(points):
    pts = sorted(set(points))
    if len(pts) <= 1:
        return pts
    def half_hull(pts):
        h = []
        for p in pts:
            while len(h) >= 2 and orient(h[-2], h[-1], p) <= 0:
                h.pop()
            h.append(p)
        return h
    lower = half_hull(pts)
    upper = half_hull(pts[::-1])
    return lower[:-1] + upper[:-1]

6.3 Monotone chain (Andrew) — O(n log n)

L'implémentation ci-dessus EST le monotone chain : tri lexicographique, puis construction de la moitié inférieure et supérieure. C'est la méthode recommandée (robuste, simple).

6.4 Jarvis march (gift wrapping) — O(n·h)

1. Commencer par le point le plus à gauche.
2. À chaque étape, choisir le point qui fait l'angle minimal (orientation max).
3. Répéter jusqu'au retour.

def convex_hull_jarvis(points):
    pts = list(set(points))
    if len(pts) <= 1:
        return pts
    hull = []
    p = min(pts)                      # point le plus à gauche
    while True:
        hull.append(p)
        q = pts[0]
        for r in pts[1:]:
            if r == p: continue
            o = orient(p, q, r)
            if o > 0 or (o == 0 and dist2(p, r) > dist2(p, q)):
                q = r
        p = q
        if p == hull[0]:
            break
    return hull

• Complexité : O(n·h) où h = nombre de points du hull. Utile si h petit.

6.5 Comparaison

---

Partie VII — Closest pair (divide & conquer)

Déjà traité en détail au chapitre 12 :

1. Trier par x, diviser en 2.
2. δ = min(closest gauche, closest droite).
3. Bande de largeur 2δ : chaque point compare avec ≤ 7 voisins triés par y.

import math

def closest_pair(points):
    pts = sorted(points)
    return _cp(pts)

def _cp(pts):
    if len(pts) <= 3:
        return min((dist2(pts[i], pts[j]) for i in range(len(pts))
                    for j in range(i + 1, len(pts))), default=float("inf"))
    mid = len(pts) // 2
    d = min(_cp(pts[:mid]), _cp(pts[mid:]))
    strip = [p for p in pts if (p[0] - pts[mid][0]) ** 2 < d]
    strip.sort(key=lambda p: p[1])
    for i in range(len(strip)):
        for j in range(i + 1, len(strip)):
            if (strip[j][1] - strip[i][1]) ** 2 >= d:
                break
            d = min(d, dist2(strip[i], strip[j]))
    return d

O(n log n).

---

Partie VIII — Line sweep

8.1 Principe

Un balayage (sweep line) fait passer une ligne verticale de gauche à droite et maintient une structure de données des événements. Événements triés → traitement en O((n+k) log n).

8.2 Exemple : intersections de segments

1. Événements : début (start), fin (end), croisement.
2. Ordre vertical des segments maintenu dans un arbre équilibré.
3. Quand deux segments se croisent → échanger leur ordre → détection de l'intersection.

8.3 Interval intersections / Union d'intervalles

def union_length(intervals):
    intervals.sort()
    total, cur_start, cur_end = 0, None, None
    for s, e in intervals:
        if cur_end is None or s > cur_end:
            if cur_end is not None:
                total += cur_end - cur_start
            cur_start, cur_end = s, e
        else:
            cur_end = max(cur_end, e)
    if cur_end is not None:
        total += cur_end - cur_start
    return total

8.4 Autres applications

• Point in polygon par ray casting (sweep).
• Skyline problem (silhouette de bâtiments).
• Plus grand rectangle dans un histogramme (stack + sweep).

---

Partie IX — Transformations

9.1 Translation

p' = p + t

9.2 Rotation d'angle θ (autour de l'origine)

x' = x·cos(θ) − y·sin(θ)
y' = x·sin(θ) + y·cos(θ)

9.3 Scaling

x' = x·sx,  y' = y·sy

9.4 Matrices homogènes

En 2D, on utilise des matrices 3×3 homogènes pour composer translation + rotation + scaling en une seule multiplication :

[ x' ]   [ a b tx ] [ x ]
[ y' ] = [ c d ty ] [ y ]
[ 1  ]   [ 0 0 1  ] [ 1 ]

La composition = produit de matrices → transformations complexes en O(1) par point.

---

Partie X — Applications

---

Récapitulatif des complexités

---

Check-list de maîtrise

• Je sais calculer un produit vectoriel et en interpréter le signe.
• Je sais déterminer l'orientation de 3 points (CCW/CW/collinéaires).
• Je sais tester l'intersection de segments (cas dégénérés inclus).
• Je sais calculer l'aire d'un polygone (shoelace).
• Je sais implémenter point-in-polygon.
• Je sais construire un convex hull (monotone chain + Jarvis).
• Je sais résoudre closest pair en O(n log n).
• Je sais appliquer une transformation géométrique.