Cercle d’Andres et application aux rotations discrètes.

1^{o}) Les méthodes classiques de tracé de cercles discrets. Le problème posé est d’approcher le tracé d’un cercle de centre à coordonnées entières et de rayon entier par celui d’une série de points à coordonnées entières (les pixels de l’écran). On souhaite que ce tracé soit le plus fidèle possible (en un sens à préciser) et qu’il donne une courbe au moins 8-connexe, c’est-à-dire que pour chaque point du cercle discret, deux de ses 8 voisins au moins appartiennent également au cercle. Rappelons ici brièvement les méthodes par DDA (Differential Digital Analyser) et de Brésenham. On s’intéresse (à une translation près) au tracé d’un cercle de centre O . Soit r son rayon (entier). Les deux méthodes procèdent par tracé d’un octant (celui correspondant à x ≥ 0 , y ≥ 0 , y ≥ x) puis par symétries, et sont incrémentales à partir du point (0, r) . En effet, dans le premier octant, la tangente au cercle est de pente en valeur absolue ≤ 1 , si bien qu’on progresse de manière connexe le long du cercle en faisant progresser l’abscisse de 1 . La question revient donc à déterminer à quelles abscisses il faut faire décroître y de 1 .

a) Méthode par analyse différentielle discrète (DDA) ; Ici on trace en respectant au mieux l’équation différentielle

x.dx + y.dy = 0 .

En fait, dx = 1 est constant, et les |dy| s’accumulent jusqu’à atteindre et dépasser, au choix, 1 ou 1/2 . Cela donne plusieurs algorithmes possibles, sur le modèle suivant :

let cercle dda r x0 y0 =

let x = ref 0 and y = ref r

(* 1 *)  and d= ref (-r)

in

while ( !x)<=( !y) do

point ( !x+x0) ( !y+y0) ; point ( !y+x0) ( !x+y0) ; x:=( !x)+1 ;

(* 2 *)  d:=( !d)+2*( !x) ;
(* 3 *)  if ( !d)>=0 then

begin

y:=( !y)-1 ;

(* 4 *)  d:=( !d)-2*( !y) ;

end ;

done ; ;

La variable d initialisée ligne (* 1 *) accumule les 2 ∗ x . Quand elle dépasse la valeur 0 , c’est que |dy| dépasse la valeur 1/2 ; on fait alors comme si dy valait −1 , c’est-à-dire qu’on décrémente y de 1 et d de −2 ∗ y (ligne (* 4 *)) de manière à respecter l’inéquation d = ∑(2x.dx + 2y.dy) ≤ 0 . Les variantes concernent la ligne (* 3 *) où on peut mettre ( !d)>0 à la place de ( !d)>=0 , et les lignes (* 2 *), (* 4 *) qu’on peut remplacer respectivement par

(* 2 *)  ... d:=( !d)+( !x) ;
(* 4 *) ... d:=( !d)-( !y) ;

Les quatre quarts de cercle obtenus sont montrés en figure 1.

b) La Méthode de Brésenham. La méthode précédente offre peu de justifications théoriques, de plus, quelle variante choisir ? Son intérêt est qu’elle peut se généraliser au tracé d’autres courbes algébriques, des coniques par exemple.

Avec la méthode de Brésenham, on a une justification théorique : on choisit à chaque étape, entre les points (x, y) et (x, y − 1) celui qui est le plus proche du point théorique (x, \sqrt{r^{2} − x^{2}}) . Or,

y − \sqrt{r^{2} − x^{2}} ≤ \sqrt{r^{2} − x^{2}} − (y − 1) ⇔ 1 + 4(x^{2} + y^{2} − y − r^{2}) ≤ 0 .

Le choix se fait donc au moyen de la variable de décision d = 1 + 4(x^{2} + y^{2} − y − r^{2}) . Cette variable se calcule de façon incrémentale à partir de sa valeur initiale 1 − 4r , d’où l’algorithme suivant :

let cercle bres r x0 y0 =

let x = ref 0 and y = ref r and d= ref (1-4*r) in

(* 1 *)  while ( !x)<=( !y) do
(* a *)  point ( !x+x0) ( !y+y0) ;
(* b *)  point ( !y+x0) ( !x+y0) ;

x:=( !x)+1 ;

d:=( !d)+8*( !x)-4 ;

(* 2 *)  if ( !d)>0 then

begin d:=( !d)-8*( !y)+8 ; y:=( !y)-1 ;

end ;

done ; ;

A la ligne (* 2 *) on peut remplacer d > 0 par d ≥ 0 ce qui change le choix et fait préférer y−1 à y en cas d’équidistance des deux points (x, y) et (x, y −1) au point théorique. Cela n’a aucune influence si le centre est à coordonnées entières, mais peut en avoir si le centre est à coordonnées demi-entières. On montre en Figure 1. différents cercles obtenus par DDA et par la méthode de Brésenham. Les cercles de Brésenham de droite sont obtenus par une modification de l’algorithme expliquée plus loin.

Remarques :

- On voit (cas des cercles de Brésenham de rayons 4 et 11 ) que, dans le programme précédent, le test ligne (* 1 *) n’est pas le meilleur possible, car il met parfois un point (x, x) sur le cercle alors que le point (x − 1, x) y est déjà, ce qui rajoute une “pointe” peu esthétique. On corrige ce défaut en permutant, dans la boucle while , l’affichage (lignes (* a*) et (* b *)) avec les lignes suivantes de calcul des nouveaux x, y, d , en rajoutant cet affichage avant la boucle while , et en remplaçcant le test de la ligne (* 1 *) par (∗1′∗) while(!x) + 1 ≤ (!y) − 1 .

- Un logiciel comme Paint Shop Pro par exemple, ou wcaml74 , utilise Brésenham comme routine de tracé de cercles.

2^{o}) Le cercle d’Andres. Le défaut des cercles précédents est qu’ils oublient des points : il y a environ 10% de trous, c’est-à-dire de points n’appartenant à aucun cercle discret de centre entier donné. En revanche, tout point (x, y) est dans une couronne r − \frac{1}{2} ≤ \sqrt{x^{2} + y^{2}} < r + \frac{1}{2} lorsque r décrit ℕ . D’où l’idée de la définition suivante

a) Définition : Le cercle discret d’Andres de centre (x_{0}, y_{0}) et de rayon r est l’ensemble des points (x, y) ∈ ℤ^{2} vérifiant les inéquations diophantiennes :

(A) (r − \frac{1}{2})^{2} ≤ (x − x_{0})^{2} + (y − y_{0})^{2} < (r + \frac{1}{2})^{2} .

Deux remarques à propos de cette définition :

- tout d’abord, on peut évidemment remplacer les intervalles [ r − \frac{1}{2}, r + \frac{1}{2} [ par tout ensemble d’intervalles de longueur 1 contigüs, en particulier les [r, r + 1[ (on obtient ainsi les cercles dits “naïfs”) ;

- le centre (x_{0}, y_{0}) d’un tel cercle est à coordonnées entières ou demi-entières. On verra en effet que les propriétés sont fondamentalement les mêmes pour les deux types de centres, de plus prendre un centre de coordonnées (1/2, 1/2) par exemple donne une plus grande symétrie au cercle et simplifie les algorithmes (par exemple, le nombre total de points d’un cercle est toujours exactement 4 fois celui d’un quadrant).

b) Propriétés des cercles d’Andres. Proposition : Le cercle d’Andres est 8-connexe et les points du cercle peuvent être tous parcourus dans un ordre canonique, donné par un algorithme de type incrémental. Pour démontrer cette proposition, plaçons-nous en un point (x, y) ∈ ℕ^{2} du cercle, et cherchons à déterminer ses voisins éventuels à droite et en-dessous, suivant la figure ci-dessous :

Lemmes :

1. Si (x, y) vérifie (A) , ni (x + 1, y + 1) ni (x − 1, y − 1) ne vérifient (A) (i.e. les cas (d) et (e) sont exclus).

2. Si (x+1, y) vérifie (A) , (x, y−1) ne le vérifie pas, et réciproquement (i.e. les cas (a) et (b) s’excluent mutuellement).

3. Pour y ≥ 1 , si ni (x + 1, y) ni (x, y − 1) ne vérifient (A) , alors (x + 1, y − 1) vérifie (A) (i.e. les cas non (a) et non (b) impliquent le cas (c) ).

Démonstrations : Les démonstrations de ces trois lemmes se font dans le cas x_{0} = y_{0} = 0 et reposent toutes sur le fait qu’on a alors

x + y ≥ \sqrt{x^{2} + y^{2}} ≥ (r − \frac{1}{2}) .

Pour le lemme 1 ,on voit facilement que

x^{2} + y^{2} ≥ (r − \frac{1}{2})^{2} ⟹ (x + 1)^{2} + (y + 1)^{2} > (r + \frac{1}{2})^{2} ,

d’où on déduit que si (x, y) vérifie (A), il est impossible que (x − 1, y − 1) vérifie (A) aussi. Démonstration du lemme 2 : Supposons (absurde) à la fois (a) et (b) . On aurait alors :

(r − \frac{1}{2})^{2} ≤ x^{2} + (y − 1)^{2} et (x + 1)^{2} + y^{2} < (r + \frac{1}{2})^{2} ,

d’où, successivement :

r^{2} − r + \frac{1}{4} < r^{2} + r + \frac{1}{4} − 2y − 2x

⟹x + y < r
⟹\sqrt{x^{2} + y^{2}} ≤ x + y ≤ r − 1  (on est avec des entiers).

Or ceci contredit le fait que (x, y) vérifie (A) , cqfd. Enfin, pour le lemme 3 : L’hypothèse est que (x + 1)^{2} + y^{2} ≥ (r + \frac{1}{2}) et x^{2} + (y − 1)^{2} < (r − \frac{1}{2})^{2} . On en déduit :

\begin{aligned}(r + \frac{1}{2})^{2} − 2y + 1 ≤ (x + 1)^{2} + (y − 1)^{2} & < (r − \frac{1}{2})^{2} + 2x + 1 \\ ⇔ (r − \frac{1}{2})^{2} + 2(r − y) ≤ (x + 1)^{2} + (y − 1)^{2} & < (r + \frac{1}{2})^{2} + 2(x − r) + 1 . \end{aligned}

Or, les équations (A) impliquent que x ≤ r nécessairement, mais si x = r , on voit de plus que (x, y − 1) vérifie (A) aussi, or c’est exclu par hypothèse. On a donc x < r et comme de même r − y ≥ 0 , on en déduit que

(r − \frac{1}{2})^{2} < (x + 1)^{2} + (y − 1)^{2} ≤ (r + \frac{1}{2})^{2} , cqfd.

On déduit de ces trois lemmes que l’algorithme initialisant (x, y) à (x_{0}, r + y_{0}) et faisant progresser le point (x, y) du cercle d’Andres suivant les cas (exclusifs) (a) ou (b) quand c’est possible, sinon suivant le cas (c) , et ceci tant que y reste ≥ 0 , fait parcourir à (x, y) tous les points du cercle d’Andres situés dans le premier quadrant. Ceci se démontre par invaraint de boucle en utilisant la relation d’ordre partiel sur les couples

(x_{1}, y_{1}) ≺ (x_{2}, y_{2}) ⇔ x_{1} < x_{2} ou (x_{1} = x_{2} et y_{1} ≥ y_{2}) .

En effet, raisonnons dans le cas du centre (0, 0) . Pour x = 0 , (0, r) est le seul point du cercle situé dans le premier quadrant et ≺ (0, r) (clair sur (A) où x = 0 ). Supposons que l’algorithme ait tracé, jusqu’au point (x, y) , tous les points ≺ (x, y) appartenant au cercle d’Andres et situés dans le premier quadrant. Le point (x′, y′) suivant choisi par l’algorithme possède, vu les 3 lemmes, aussi cette propriété : si c’est le point choisi par (b) c’est évident puisque les points précédents (x′, y′) par ≺ sont les prédécesseurs de (x, y) et (x, y) lui-même ; si c’est le point choisi par (a) , comme le point correspondant à (e) n’est pas sur le cercle, et donc tous ceux qui sont au-dessus a fortiori non plus, de même que le point choisi par (b) (lemme 2) et a fortiori tous les points en-dessous, c’est encore vrai ; enfin, si c’est le point choisi par (c) , comme ne sont pas sur le cercle les points correspondant à (a) , et a fortiori tous ceux qui sont au-dessus, ni ceux correspondant à (b) et a fortiori tous ceux qui sont en-dessous, c’est toujours vrai. Or, l’algorithme finit nécessairement par atteindre le cas y = 0 , x = r , puisque le choix (a) ne peut être effectué au plus que r fois et que pour les autres choix, possibles tant que y ≥ 1 , y diminue de 1 . Ceci achève de démontrer notre proposition.

Donnons pour terminer l’algorithme de parcours du cercle d’Andres en CAML . La condition pour le choix (a) s’écrit :

(x + 1)^{2} + y2 < (r + \frac{1}{2})^{2} ⇔ r^{2} + r − x^{2} − y^{2} − 1 ≥ 2x

(on tient compte du 1/4 en autorisant l’égalité, puisqu’on est avec des entiers), et celle pour le choix (b) est

x^{2} + (y − 1)^{2} ≥ (r − \frac{1}{2})^{2} ⇔ r^{2} + r − x^{2} − y^{2} − 1 < 2(r − y) ,

(ici on peut ne plus tenir compte du 1/4 ). Comme le cas (c) est choisi automatiquement quand les cas (a) et (b) font défaut, on voit que tout peut se faire à l’aide de la seule variable de contrôle d = r^{2} + r − x^{2} − y^{2} − 1 . De plus, celle-ci se calcule par incrémentation au fur et à mesure à partir de sa valeur initiale d = r − 1 . Algorithme : (* parcours du quart de cercle d’Andres de centre (x_{0}, y_{0}) ∈ ℤ^{2} *)

let cercle a r x0 y0 =

let x = ref 0 and y = ref r and d= ref (r-1) in

while 0<=( !y) do

point ( !x+x0) ( !y+y0) ;

(* (a) *)  if ( !d)>=2*x

then (d:=( !d)-2*( !x)-1 ; x:=( !x)+1)

else

(* (b) *)  if ( !d)<2*(r-y)

then (d:=( !d)+2*( !y)-1 ; y:=( !y)-1)

(* (c) *)  else

(d:=( !d)+2*(( !y)-( !x)-1) ; x:=( !x)+1 ; y:=( !y)-1)

done ; ;

Remarque : On ne peut prolonger cette approche de manière naturelle pour les ellipses. En effet, soit a et b les deux demi-axes (entiers) d’une ellipse, et posons a = \frac{p}{q}.b , où \frac{p}{q} ≥ 1 est une fraction irréductible. L’équation de l’ellipse se ramène à :

q^{2}x^{2} + p^{2}y^{2} = q^{2}a^{2} = p^{2}b^{2} .

Si, par analogie avec les cercles on pose comme équations de l’ellipse discrète

q^{2}(a − \frac{1}{2})^{2} ≤ q^{2}x^{2} + p^{2}y^{2} < q^{2}(a + \frac{1}{2})^{2} (1) ou p^{2}(b − \frac{1}{2})^{2} ≤ q^{2}x^{2} + p^{2}y^{2} < p^{2}(b + \frac{1}{2})^{2} (2) ,

alors ou bien l’ellipse est non connexe (équation (1) ) ou bien il y a trop de points sur l’ellipse si bien que la propriété d’ordre naturel sur les points n’est plus satisfaite (équation (2) ). Par exemple, pour a = 30, b = 20 on trouvera les ensembles de points suivants :

Equation (1)  Equation (2)

c) Cercles de centre demi-entier. On s’intéressera ici aux cercles de centre (1/2, 1/2) , les autres cas étant similaires. L’équation est donc :

(A′) (r − \frac{1}{2})^{2} ≤ (x − \frac{1}{2})^{2} + (y − \frac{1}{2})^{2} < (r + \frac{1}{2})^{2} .

On laisse au lecteur la démonstration d’une proposition strictement analogue au cas précédent, justifiant le tracé de tous les points entiers vérifiant (A′) au moyen d’un algorithme incrémental. Soit (x, y) un point vérifiant (A′) : le choix du point suivant (x + 1, y) (a) se fait si

(x + 1 − \frac{1}{2})^{2} + (y − \frac{1}{2})^{2} < (r + \frac{1}{2})^{2} ⇔ r^{2} + r − x^{2} − y^{2} + x + y > 2x ,

le choix du point suivant (x, y − 1) (b) se fait si

(x − \frac{1}{2})^{2} + (y − 1 − \frac{1}{2})^{2} ≥ (r − \frac{1}{2})^{2} ⇔ r^{2} + r − x^{2} − y^{2} + x + y ≤ 2(r + 1 − y) ,

et sinon on choisit le point suivant (x + 1, y − 1) (c) . On remarque à ce moment que la variable de décision d = r^{2} + r − x^{2} − y^{2} + x + y est toujours paire et doit se comparer à des nombres pairs, donc on peut simplifier les tests par 2 . Finalement, l’algorithme est le suivant : Algorithme : (* parcours du quart de cercle d’Andres de centre (x_{0} + 1/2, y_{0} + 1/2) , (x_{0}, y_{0}) ∈ ℤ^{2} *)

let cercle a r x0 y0 =

let x = ref 0 and y = ref r and d= ref r in

while 0 <=( !y) do

point ( !x+x0) ( !y+y0) ;

(* (a) *)  if ( !d)>x

then (d:=( !d)-( !x) ; x:=( !x)+1)

else

(* (b) *)  if ( !d)<=r+1-y

then (d:=( !d)+( !y)-1 ; y:=( !y)-1)

(* (c) *)  else

(d:=( !d)+( !y)-( !x)-1 ; x:=( !x)+1 ; y:=( !y)-1)

done ; ;

3^{o}) Applications aux rotations discrètes.

a) Généralités. Les propriétés précédentes du cercle d’Andres permettent de définir des coordonnées polaires discrètes. En effet, tout point appartient à un cercle de centre O(0, 0) et un seul, et sur ce cercle, on dipose d’un sens de parcours permettant de numéroter un à un tous les pixels. Dans cette optique, on définit la longueur d’un cercle comme le nombre de points le composant, et l’angle polaire comme le quotient entier entre le nombre de points à parcourir depuis le point sur l’axe polaire et la longueur du cercle. Mais des difficultés aparaissent : d’une part, la longueur n’est pas donnée par une fonction simple du rayon, mais doit être calculée algorithmiquement. On s’aperçoit à cette occasion que le quotient de la longueur sur le rayon oscille entre 6.07 et 6.50 pour des rayons de l’ordre de 200 . Mais de plus, la progression d’un point le long d’un cercle ne se fait pas à la même vitesse angulaire d’un cercle à l’autre, et pire, ne se fait pas même à vitesse angulaire constante sur un même cercle. On peut s’en douter à l’aspect tantôt lisse, tantôt “frippé” des arcs de cercles, mais aussi en représentant l’ensemble des points ayant même angle polaire. Sur la planche 1 (fichier image “arcs3030.pcx”) on a représenté des points d’angle polaire correspondant à 30^{o} , puis les points d’angles polaires décalés de 30^{o} par rapport à ceux-ci (l’axe polaire étant l’axe vertical). On voit que les droites discrètes obtenues ne sont même pas 8-connexes !

Quant à l’additivité des angles, elle est mise en défaut uniquement à cause de l’approximation dans le quotient entier par la longueur. Parfois, cela se passe bien, comme en témoigne la planche 2 (“arcs3060.pcx”) où on retrouve l’axe horizontal après une rotation de 30^{o} + 60^{o} . b) Implémentation d’une rotation discrète.

Malgré ces inconvénients (pré)visibles, nous avons implémenté voir une rotation discrète d’images en CAML . Pour cela, nous nous sommes servis de l’algorithme de parcours précédent pour des cercles de centre à coordonnées demientières, et nous en avons fait une fonction

pix suiv : int −> int −> int ref −> int −> int ∗ int ,

où pix suiv x y d r donne le pixel suivant (x, y) sur le cercle de centre (1/2, 1/2) de rayon r tout en modifiant la variable de décision d correspondante. La rotation elle-même se fait alors par parcours simultanés et décalés d’un angle polaire constant de 2 cercles de même rayon croissant, l’un dans l’image de départ, l’autre dans l’image à obtenir par rotation, et en recopiant les couleurs d’une image à l’autre.

Le résultat est visible aux planches 3 et 4 . En planche 3, on a fait subir à l’image une rotation de 135^{o} , et en planche 4 on a rajouté une rotation de 225^{o} : là encore on remarque une non-additivité des arcs discrets.

Certes, l’image est déformée, mais tout de même reconnaissable, et de plus celle-ci est obtenue sans perte d’informations. Avec un logiciel comme Paint Shop Pro, si on tourne une image de 135^{o} puis de 225^{o} on obtient une image qui n’est plus identique à l’image de départ, mais ayant subie un léger lissage (cf comparaison à l’écran des images “vis p0.pcx” et “vis p360.pcx”).

⋇ ⋇ ⋇