L'algorithme de cercle de Bresenham

Un autre cercle

Dans un précédent article j'ai parlé d'un moyen de dessiner un cercle en utilisant des fonctions trigonométriques:

		 $x = x_{C} + r * cos θ$ 

		 $y = y_{C} + r * sin θ$

Où

(x_{C}, y_{C})

sont les coordonnées du centre et r est le rayon.
Cette fois on va utiliser une autre formule pour décrire un cercle.
Rappelez-vous que par définition un cercle est l'ensemble des points qui sont à une distance donnée du centre.
Cette distance est le rayon.
En utilisant le théorème de Pythagore, on peut écrire:

		 ${(x - x_{C})}^{2} + {(y - y_{C})}^{2} = r^{2}$

Ou en prenant la racine carrée:

		 $\sqrt{{(x - x_{C})}^{2} + {(y - y_{C})}^{2}} = r$

Maintenant considérons que notre cercle est centré en (0, 0).
Imaginez qu'on est en haut du cercle, où (x, y) = (0, r).
On va augmenter x vers la droite.

Le point va avoir tendance à s'éloigner du centre.
Alors à chaque coordonnée x, on va calculer

\sqrt{x^{2} + y^{2}}

.
Et si c'est plus grand que r, on va descendre d'un pixel.
Ce sont les flèches rouges sur cette image:

On s'arrêtera quand x = y.
Quand on affiche notre pixel on n'a plus qu'a ajouter les coordonnées du centre, et le code ressemblera à ça:

		int cX = SCREEN_WIDTH  / 2;
		int cY = SCREEN_HEIGHT / 2;
		int r = 150;

		int x = 0;
		int y = r;

		while (x <= y)
		{
			gfx.setPixel(cX + x, cY + y, Color(255, 255, 255));

			x++;

			if (sqrt(x * x + y * y) > r)
				y--;
		}
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows

Et voilà le résultat.
remarquez qu'on a dessiné une partie en bas du cercle, parce qu'avec les coordonnées écran, l'axe Y augmente vers
le bas.

Maintenant on ne dessine qu'un huitième du cercle: la partie rouge dans l'image ci-dessous.
Mais on peut tirer avantage des symétries du cercle pour dessiner les autres parties en combinant simplement les
coordonnées:

Voilà ce qu'on va changer dans le code:

		while (x <= y)
		{
			gfx.setPixel(cX + x, cY + y, Color(255, 255, 255));
			gfx.setPixel(cX + x, cY - y, Color(255, 255, 255));
			gfx.setPixel(cX - x, cY + y, Color(255, 255, 255));
			gfx.setPixel(cX - x, cY - y, Color(255, 255, 255));
			gfx.setPixel(cX + y, cY + x, Color(255, 255, 255));
			gfx.setPixel(cX + y, cY - x, Color(255, 255, 255));
			gfx.setPixel(cX - y, cY + x, Color(255, 255, 255));
			gfx.setPixel(cX - y, cY - x, Color(255, 255, 255));

			x++;

			if (sqrt(x * x + y * y) > r)
				y--;
		}
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows

Et maintenant on peut voir le cercle complet.

Chronométrer

On a vu l'idée de base de la façon de tracer un cercle.
Mais comme c'était le cas pour les lignes, l'algorithme de Bresenham est plus à propos d'optimiser cette idée de
base.
Alors pour mieux voir les effets de cette optimisation, on va chronométrer cet algorithme.
D'abord on va le mettre dans une fonction:

		void circle(int centerX, int centerY, int radius, Color c)
		{
			int x = 0;
			int y = radius;

			while (x <= y)
			{
				gfx.setPixel(centerX + x, centerY + y, c);
				gfx.setPixel(centerX + x, centerY - y, c);
				gfx.setPixel(centerX - x, centerY + y, c);
				gfx.setPixel(centerX - x, centerY - y, c);
				gfx.setPixel(centerX + y, centerY + x, c);
				gfx.setPixel(centerX + y, centerY - x, c);
				gfx.setPixel(centerX - y, centerY + x, c);
				gfx.setPixel(centerX - y, centerY - x, c);

				x++;

				if (sqrt(x * x + y * y) > radius)
					y--;
			}
		}

Ensuite, pour plus de précision, on va chronométrer 10000 appels à cette fonction.

		// time 10000 circles
		sys.StartPerfCounter();

		for (int i = 0; i < 10000; i++)
			circle(SCREEN_WIDTH  / 2, SCREEN_HEIGHT / 2, 150, Color(255, 255, 255));

		printf("%f\n", sys.StopPerfCounter());
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows

Ce code prend environ 33 ms sur mon ordinateur.

Affiner les coordonnées

Dans les programmes précédents on a implicitement défini que les coordonnées entières se trouvaient à la limite
entre 2 pixels.
Ca cause quelques imperfections. Si vous regardez le pixel le plus en haut par exemple, il semble qu'il soit seul
sur sa ligne. Tous les autres pixels sont au moins une ligne en-dessous.

L'algorithme classique de Bresenham par contre considère que les coordonnées représentent le centre d'un pixel.
Disons que nous sommes au point (x, y). Le prochain point que nous allons dessiner est soit celui à droite (x+1, y),
soit celui qui se trouve en-dessous (x+1, y-1)

On va seulement tester un point entre ces deux là, aux coordonnées (x+1, y-0.5).

Comme avant, on va tester si:

		 $\sqrt{x^{2} + y^{2}} > r$

Mais on va modifier un peu ça.
D'abord, la racine carrée est une opération coûteuse. alors on peut aussi bien tester le carré de ces valeurs:

		 $x^{2} + y^{2} > r^{2}$

Ensuite, pour se rapprocher de l'algorithme original, on va soustraire r² des deux cotés:

		 $x^{2} + y^{2} - r^{2} > 0$

Le code à la fin de la boucle ressemble maintenant à ca:

		float X = x + 1;
		float Y = y - 0.5f;
		float m = X * X + Y * Y - radius * radius;

		x++;

		if (m > 0)
			y--;

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Ce code n'est pas vraiment plus rapide, mais il produit un cercle un peu plus rond.

Utiliser des entier

Les calculs sur des entiers sont plus rapides que sur des nombres flottants.
Mais les équations qu'on utilise ne sont pas vraiment adaptées aux entiers.
Alors retravaillons-les. On a:

		 $X = x + 1$ 

		 $Y = y - 0.5$ 

		 $m = X^{2} + Y^{2} - r^{2}$

Si on remplace X et Y dans la 3ième équation, on obtient:

		 $m = {(x + 1)}^{2} + {(y - 0.5)}^{2} - r^{2}$

Maintenant si on développe ça

		 $m = x^{2} + 2 x + 1 + y^{2} - y + 0.25 - r^{2}$ 

		 $m = x^{2} + y^{2} + 2 x - y - r^{2} + 1.25$

Maintenant multiplions tout par 4:

		 $4 m = 4 (x^{2} + y^{2} + 2 x - y - r^{2}) + 5$

Appelons cette valeur "M".
Comme x, y et r sont des entiers, M est aussi un entier.
Et comme on teste seulement si m est plus grand que 0, le résultat serait le même si on testait 4m > 0 ou M > 0.
Donc notre code va maintenant ressembler à ça:

		int m = 4 * (x*x + y*y + 2*x - y - radius*radius) + 5;

		x++;

		if (m > 0)
			y--;

		Télécharger le code source
		Télécharger l'exécutable pour Windows

Cette version prend environ 10 ms.

Incrémentation

Comme pour l'algorithme de ligne, le principe central de celui de cercle est d'incrémenter une valeur d'erreur pour
chaque pixel qu'on dessine.
Ici "m" peut être considéré comme notre valeur d'erreur.
Pour l'instant, elle est recalculée à chaque étape à partir des valeurs de x, y et r.
Mais on peut imaginer de l'initialiser avec une valeur de départ, puis de lui ajouter une certaine valeur à chaque
fois qu'on incrémente x, et une autre valeur à chaque fois qu'on décrémente y.

Cette explication est peut-être confuse. Alors voyons ça dans le code.
Pour le moment notre boucle ressemble à ça:

		int x = 0;
		int y = radius;

		while (x <= y)
		{
			[draw pixels]

			int m = 4 * (x*x + y*y + 2*x - y - radius*radius) + 5;

			x++;

			if (m > 0)
				y--;
		}

Et on veut la modifier pour qu'elle ressemble à ça:

		int x = 0;
		int y = radius;
		int m = m_start;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				m += m_inc_y;
				y--;
			}

			m += m_inc_x;
			x++;
		}

Et tout ce qu'on a à faire c'est de trouver les valeurs de m_start, m_inc_x et m_inc_y.

Alors essayons de trouver m_inc_x.
Au point courant on a:

		 $m = 4 (x^{2} + y^{2} + 2 x - y - r^{2}) + 5$

Si on incrémente x de 1, on a:

		 $m' = 4 ({(x + 1)}^{2} + y^{2} + 2 (x + 1) - y - r^{2}) + 5$

Travaillons ça:

		 $m' = 4 (x^{2} + 2 x + 1 + y^{2} + 2 x + 2 - y - r^{2}) + 5$ 

		 $m' = 4 (x^{2} + y^{2} + 2 x - y - r^{2}) + 5 + 4 (2 x + 1 + 2)$ 

		 $m' = m + 4 (2 x + 3)$ 

		 $m' = m + 8 x + 12$

Donc

		 $m_inc_x = 8 x + 12$

Et notre code devient:

		int x = 0;
		int y = radius;
		int m = m_start;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				m += m_inc_y;
				y--;
			}

			m += 8 * x + 12;
			x++;
		}

Alternativement on peut augmenter "m" après avoir incrémenté x. Dans ce cas on a:

		int x = 0;
		int y = radius;
		int m = m_start;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				m += m_inc_y;
				y--;
			}

			x++;
			m += 8 * x + 4;
		}

Maintenant pour m_inc_y. Quand on décrémente y de 1 on a:

		 $m' = 4 (x^{2} + {(y - 1)}^{2} + 2 x - (y - 1) - r^{2}) + 5$ 

		 $m' = 4 (x^{2} + y^{2} - 2 y + 1 + 2 x - y + 1 - r^{2}) + 5$ 

		 $m' = 4 (x^{2} + y^{2} + 2 x - y - r^{2}) + 5 + 4 (- 2 y + 1 + 1)$ 

		 $m' = m + 4 (- 2 y + 2)$ 

		 $m' = m - 8 y + 8$ 

		 $m_inc_y = - 8 y + 8$

Ca mène à ce code:

		int x = 0;
		int y = radius;
		int m = m_start;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				m += -8 * y + 8;
				y--;
			}

			x++;
			m += 8 * x + 4;
		}

Ou si on augmente m après avoir décrémenté y:

		int x = 0;
		int y = radius;
		int m = m_start;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				y--;
				m -= 8 * y;
			}

			x++;
			m += 8 * x + 4;
		}

Maintenant pour m_start, pour rappel, on a:

		 $m = 4 (x^{2} + y^{2} + 2 x - y - r^{2}) + 5$

Et on est au point de coordonnées (0, rayon).
Donc si on remplace ça dans l'équation, on obtient:

		 $m = 4 (0^{2} + r^{2} + 2 * 0 - r - r^{2}) + 5$ 

		 $m = 4 (r^{2} - r - r^{2}) + 5$ 

		 $m = -4 r + 5$

Et le code final est:

		int x = 0;
		int y = radius;
		int m = 5 - 4 * radius;

		while (x <= y)
		{
			[draw pixels]

			if (m > 0)
			{
				y--;
				m -= 8 * y;
			}

			x++;

			m += 8 * x + 4;
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows

C'est fondamentalement le même niveau d'optimisation que l'algorithme historique de Bresenham's.
Mais bien que les calculs à l'intérieur de la boucle sont beaucoup moins complexes qu'avant, cette version du code
n'est pas plus rapide que les programme précédent sur un PC actuel.
C'est probablement à cause de la façon dont le compilateur optimise le code et de la façon dont les processeurs
modernes fonctionnent.

Dans le prochain article on verra comment dessiner des cercles pleins à partir de cet algorithme et comment les
optimiser.