Aller au contenu

Géométrie Discrète - CM séance 03

3. Approximation polygonale d'un contour

Objectif : sélectionner parmi les points du contours une sous-liste de sommets d'un polygone, qui approxime "au mieux" le contour.

À l'issue de l'opération, chaque point du contour sera à une distance du polygone qui sera inférieure à un seuil fixé seuil_dist.

3.1 Distance à une droite

        A                            A                  D
        +                            +------------------+ 
        |                           /|                 /  
        |d(A,(BC))                 / |h               /   
        |                         /  |               /    
  --+------------------+--       +------------------+     
    B                  C         B                  C

On veut calculer \(d(A,(BC))\).

Soit P le parallélogramme \((B,C,D,A)\). On a :

  • aire(P) = \(||BC|| \times h\), où \(h\) = \(d(A,(BC))\) ;
  • on a aussi : aire(P) = \(|\text{det}(\vec{BC}, \vec{BA})|\)

Donc \(d(A,(BC))\) = \(h\) = \(|\text{det}(\vec{BC}, \vec{BA})| \,/\, ||BC||\)

\[ = \frac{ \text{abs} \left| \begin{array}{cc} x_C-x_B & x_A-x_B \\ y_C-y_B & y_A-y_B \end{array} \right| }{ \sqrt{(x_C-x_B)^2+(y_C-y_B)^2} } \]

\(= \text{abs}\,((x_C-x_B)*(y_A-y_B)-(y_C-y_B)*(x_A-x_B))/\sqrt{(x_C-x_B)^2+(y_C-y_B)^2}\)

si B et C sont distinct.

Si B et C sont confondus, alors \(d(A,(BC))\) = \(d(A,B)\) = \(\sqrt{(x_B-x_A)^2+(y_B-y_A)^2}\).

Dans notre approximation, chaque segment du polygone représente un tronçon du contour ; donc chaque point du contour devra satisfaire le seuil de distance pour la droite définie par le segment qui représente son propre tronçon.

Autrement dit : pour chaque point du contour, on ne testera sa distance qu'à un seul segment du polygone.

3.2 Algorithme récursif avec seuil

Algorithme de Ramer-Douglas-Peucker :

  • [Ramer 1972] An iterative procedure for the polygonal approximation of plane curves. CGIP.
  • [Douglas, Peucker 1973] Algorithms for the reduction of the number of points required to represent a digitized line or its caricature. The Canadian Cartographer.

En entrée :

  • une chaîne de Freeman d'un contour cfc[]
  • les coordonnées du point de départ \(A\) du contour

Initialisation :

  • on calcule, à partir des coordonnées de \(A\) et des directions de cfc[], les coordonnées des points du contour, que l'on stocke dans des tableaux cont_x[], cont_y[], y compris pour le dernier point \(A' = A\) que l'on rajoute pour fermer le contour ;
  • on mémorise aussi pour chaque point un flag initialisé à 1 dans cont_f[] : ce flag signifie "est un sommet du polygone".

L'algorithme va consister à mettre la plupart des flags à 0, de façon à ne conserver que les sommets du polygone recherché.

Dans tout l'algorithme, on désigne les points par leur indice dans les tableaux.

Étape 1 :

  • on cherche le point (noté \(B\)) le plus éloigné de \(A\) ;
  • on lance l'appel récursif sur \([ A \ldots B ]\) et sur \([ B \ldots A']\).

Étape 2 : appel récursif pour un tronçon \([P \ldots Q]\) :

  • si le tronçon est de taille \(\leq 2\) on sort ;
  • on cherche le point \(R\) du tronçon qui est le plus éloigné de la droite \((PQ)\) ;
  • si \(d(R,(PQ)) <\) seuil_dist on met le flag "sommet" cont_f à 0 sur les points \(]P \ldots Q[\) ;
  • sinon on fait un appel récursif sur \([P \ldots R]\) et sur \([R \ldots Q]\).

Étape 3 : relaxation

  • on parcourt une seule fois la liste des sommets : si pour un triplet de sommets consécutifs \(ABC\) on a \(d(B,(AC)) <\) seuil_dist, alors on met le flag "sommet" de \(B\) à 0 dans cont_f.

4. Droites et plans discrets

4.1 Algorithme de Bresenham

[Bresenham 1965] Algorithm for computer control of a digital plotter. IBM Systems Journal.

Algorithme de tracé d'un segment discret en nombres entiers :

  • Pour une pente entre 0 et 1, c'est-à-dire dans le premier octant :
def calculer_segment (x0, y0, x1, y1) :
    dx = x1 - x0
    dy = y1 - y0
    assert 0 <= dy <= dx
    #
    k = 2*dy - dx               # variable de contrôle
    y = y0
    for x in range(x0, x1+1) :  # de x0 à x1
        print(x, y, k)          # dessiner (x,y)
        if k > 0 :
            y += 1
            k -= 2*dx
        k += 2*dy
    assert y == y1
  • Pour les autres directions : par symétrie.

Exemple : calculer_segment (1, 2, 15, 5)

Observations :

  • le segment est 8-connexe par construction ;
  • la longueur des paliers est \(\leq \lceil\frac{dx}{dy}\rceil\) (arrondi supérieur) ;
  • la longueur des paliers diffère de 1, sauf les paliers extrêmes qui sont en fait des demi-paliers.

4.2 Droites de Reveillès

Différents essais de caractérisation d'une droite discrète dans la littérature :

Définition arithmétique de [Reveillès 1991, p.43] :
Soit \(a,b,c,\omega\) des entiers tels que \(\text{pgcd}(a,b) = 1\).
Une droite discrète de paramètres \((a,b,c,\omega)\) est l'ensemble des points de \(\Bbb{Z}^2\) vérifiant : \(0 \leq ax + by + c < \omega\) .

Le paramètre \(\omega\) est l'épaisseur arithmétique de la droite.

Théorème [Reveillès 1991, pp.48-49] :

  • si \(\omega < \max(|a|,|b|)\) alors la droite n'est pas connexe ;
  • si \(\omega = \max(|a|,|b|)\) alors la droite est 8-connexe (droite de Bresenham) ;
  • si \(\max(|a|,|b|) < \omega < |a|+|b|\) alors la droite est \(*\)-connexe ;
  • si \(\omega = |a|+|b|\) alors la droite est 4-connexe ;
  • si \(\omega > |a|+|b|\) la droit est dite épaisse.
Exemple : \(0 \leq 3x - 7y + 0 < \omega\)
\(\omega = 5\) : non connexe
\(\omega = 7\) : 8-connexe
\(\omega = 9\) : \(*\)-connexe
\(\omega = 10\) : 4-connexe
\(\omega = 13\) : droite épaisse
Généralisation aux plans discrets :
Soit \(a,b,c,d,\omega\) des entiers tels que \(a \leq b \leq c\), \(c \neq 0\).
Un plan discret de paramètres \((a,b,c,d,\omega)\) est l'ensemble des points de \(\Bbb{Z}^3\) vérifiant : \(0 \leq ax + by + cz + d < \omega\)
où \(\omega\) est l'épaisseur arithmétique du plan.

Les sections du plan discret avec les plans \((Oxy)\), \((Oxz)\) ou \((Oyz)\) sont alors des droites discrètes, dont on peut étudier les paliers et la connexité.

La structure topologique d'un plan discret dépend de la présence de tunnels :
un plan discret \(0 \leq ax + by + cz + d < \omega\) possède un \(k\)-tunnel \((k = 0,1,2)\) s'il existe deux \(k\)-voisins \(E\) et \(F\), l'un en dessous et l'autre au-dessus du plan, autrement dit si
\(a\,x_E + b\,y_E + c\,z_E + d < 0\)   et   \(a\,x_F + b\,y_F + c\,z_F + d >= \omega\) .

Proposition [Andrès 1993] :

  • si \(\omega < c\) alors le plan a des 2-tunnels ;
  • si \(c \leq \omega < b+c\) alors le plan a des 1-tunnels ;
  • si \(b+c \leq \omega < a+b+c\) alors le plan a des 0-tunnels ;
  • si \(a+b+c \leq \omega\) alors le plan n'a pas de tunnels.

De la présence des tunnels on déduit la connexité :

Théorème [Andrès 1993] :

  • si \(\omega = c\) alors le plan est 18-connexe ;
  • si \(c < \omega < b+c\) alors le plan est 18 ou 6-connexe ;
  • si \(b+c \leq \omega\) alors le plan est 6-connexe.

4.3 Reconnaissance de droite discrète

Algorithme arithmétique d'Isabelle Debled-Rennesson : [Debled-Renesson 1995] A linear algorithm for segmentation of digital curves. IJPRAI.

Données : une suite de \(n\) points \(P_i(x_i,y_i)\) formant un chemin 8-connexe.
Contraintes : \(P_0 = O\) ; \(x_{i+1} = x_i+1\) ; le chemin est contenu dans le 1er octant : \(0 \leqslant y_i \leqslant x_i\).
But : reconnaître un segment de droite 8-connexe en rajoutant des \(P_i\) ; déterminer son équation \(0 \leq ax -by +c < \omega = max(|a|,|b|)\) ; s'interrompre dès qu'on n'a plus un segment.

Principe : durant l'algorithme on met à jour les points d'appui inférieurs et supérieurs d'abscisses minimale et maximale \(L\), \(L'\), \(U\), \(U'\) → ceci contraint de plus en plus l'équation du segment.

# Initialisation :
a = y_1, b = 1, c = 0
U = L = (0,0)
U' = L' = (1,y_1)

# Parcours
pour i de 2 à n-1 :
    M = (x_i,y_i)
    r = a*x_i -b*y_i +c

    si 0 < r < b-1 :            # cas A : ok
        continuer
    sinon si r == 0 :           # M sur droite d'appui supérieure
        U' = M
    sinon si r == b-1 :         # M sur droite d'appui inférieure
        L' = M
    sinon si r < -1 ou r > b :  # cas B : ce n'est pas un segment
        renvoie Faux
    sinon si r == -1 :          # cas C : M juste au dessus de (UU')
        L = L' ; U' = M
        a = y_i-yU ; b = x_i-xU
        c = -a*x_i + b*y_i
    sinon si r == b :           # cas D : M juste en dessous de (LL')
        U = U' ; L' = M
        a = y_i-yL ; b = x_i-xL
        c = -a*x_i + b*y_i +b-1

# Résultat
renvoie Vrai  # c'est un segment défini par a,b,c

Exemple :

P = [(0, 0), (1, 0), (2, 1), (3, 1), (4, 1), (5, 2), (6, 3)]
init  U = (0, 0) L = (0, 0) U' = (1, 0) L' = (1, 0) a = 0 b = 1 c = 0
i = 2 M = (2, 1) -> r = -1
   C: U = (0, 0) L = (1, 0) U' = (2, 1) L' = (1, 0) a = 1 b = 2 c = 0
i = 3 M = (3, 1) -> r = 1
 b-1: U = (0, 0) L = (1, 0) U' = (2, 1) L' = (3, 1) a = 1 b = 2 c = 0
i = 4 M = (4, 1) -> r = 2
   D: U = (2, 1) L = (1, 0) U' = (2, 1) L' = (4, 1) a = 1 b = 3 c = 1
i = 5 M = (5, 2) -> r = 0
   0: U = (2, 1) L = (1, 0) U' = (5, 2) L' = (4, 1) a = 1 b = 3 c = 1
i = 6 M = (6, 3) -> r = -2
   B: False
Fin, i = 5 a = 1 b = 3 c = 1