L'algorithme de ligne de Bresenham

A propos du code

Le code dans cet article a été écrit avec Code::Blocks et la SDL 2.
Vous pouvez trouver ici un guide pour installer ces logiciels.
Bien qu'il soit basé sur la SDL, je n'utilise pas ses fonctions directement. J'ai écrit une petite bibliothèque
avec quelques fonctions basiques pour faciliter la compréhension et la portabilité dans un autre langage.
Vous pouvez en apprendre plus sur cette bibliothèque ici.

Une simple ligne

Comment votre ordinateur dessine une ligne à l'écran ?
C'est assez facile à comprendre pour une ligne verticale ou horizontale, mais qu'en est-il pour une ligne de pente
quelconque que vous pourriez tracer dans un programme de dessin ?
Observons ce diagramme.


Il montre les pixels que nous devons afficher pour représenter la ligne idéale en vert.
On peut remarquer qu'un seul pixel est allumé pour chaque coordonnée x. C'est celui qui est le plus proche de la
ligne idéale.

Pour simplifier les explication, à partir de maintenant on ne parlera que des lignes comme celle de ce diagramme,
qui respectent ces contraintes:
Nous parlerons des autres cas à la fin, une fois qu'on aura bien défini l'algorithme pour ce cas là.

Maintenant comment trouver l'équation pour dessiner cette ligne ?
Rappelons-nous une équation que tout le monde doit connaître: 
		y = a * x + b
a est facile à comprendre. c'est la pente de la droite (la quantité d'augmentation en y divisée par la quantité
d'augmentation en x).
Dans notre cas: 
		a = (y1 - y0) / (x1 - x0)
Je ne veux pas calculer b, alors retravaillons un peu notre équation. On avait: 
		y = a * x + b
En particulier, notre premier point (x0, y0) appartient à la droite, donc il suit aussi l'équation: 
		y0 = a * x0 + b
Maintenant si on soustrait cette équation à la précédente, on obtient: 
		y - y0 = a * (x - x0) + b - b
Et b disparaît. Donc après une petite réécriture, on obtient une nouvelle formule pour notre ligne: 
		y = a * (x - x0) + y0
Le code le plus simple pour dessiner notre ligne devrait donc ressembler à ça: 
		// line coordinates
		int x0 = 0;
		int y0 = 0;
		int x1 = 30;
		int y1 = 10;

		float slope = (float)(y1 - y0) / (float)(x1 - x0);

		for (int x = x0; x <= x1; ++x)
		{
			float y = slope * (x - x0) + y0;
			gfx.setPixel(x, y, Color(255, 255, 255));
		}
		gfx.render();
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
Et voilà le résultat:


Cette ligne a l'air bien mais il y a un petit quelque chose qui cloche dedans.
Si vous regardez le premier et le dernier segment horizontal, ils devraient être équilibrés.
C'est seulement un problème d'arrondi. Ca serait facile à corriger, mais laissons-ça de coté parce qu'on s'en
occupera plus tard d'une autre façon.

Donc on a une routine de dessin qui marche. Alors quel est l'intérêt de l'algorithme de Bresenham ?
Est-ce qu'il permet de faire une ligne plus belle ?
En fait, pas vraiment.

Bresenham était quelqu'un qui travaillait pour IBM. Il a écrit cet algorithme en 1962.
A l'époque les ordinateurs étaient très lents et tout devait être optimisé.

En fait, son algorithme est seulement une version très optimisée de celui qu'on vient d'écrire.
Aujourd'hui beaucoup de gens pensent que l'optimisation n'est pas utile parce que les ordinateurs son beaucoup plus
efficaces.
Alors est-ce que ça vaut vraiment la peine ? Eh bien nous le vérifierons nous-même plus tard.
Quoi qu'il en soit c'est un bon exercice pour comprendre l'optimisation et comment un processeur fonctionne.
Alors regardons ça pas à pas.

Comment l'accélérer ?

Avant de répondre à cette question, demandons-nous ce qui ralentit ce code.
Regardons encore le code: 
		// line coordinates
		int x0 = 0;
		int y0 = 0;
		int x1 = 30;
		int y1 = 10;

		float slope = (float)(y1 - y0) / (float)(x1 - x0);

		for (int x = x0; x <= x1; ++x)
		{
			float y = slope * (x - x0) + y0;
			gfx.setPixel(x, y, Color(255, 255, 255));
		}
Remarque: j'ai omis la fonction gfx.render() parce qu'elle ne fait pas partie de l'algorithme lui-même, elle est
seulement appelée à la fin pour afficher le résultat à l'écran.

Alors, que voyons-nous dans ce code ? Il y a:
Voici quelques règles importantes que l'on peut donner à propos de tous ces éléments.

Quand vous avez une boucle, tout ce qui est à l'intérieur est exécuté plus de fois que ce qui est à l'extérieur.
Dans certains programmes ça peut être vraiment beaucoup plus de fois...
Dans notre code, la boucle est exécutée 21 fois. C'est la raison pour laquelle on calcule la pente à l'extérieur.
Parce que, comme c'est une valeur constante, il n'y a pas besoin de la calculer 20 fois...

Les entiers sont plus rapides que les flottants. Même si de nos jours tous les processeurs ont un coprocesseur
arithmétique intégré, les calculs sur des entiers sont toujours plus rapides que sur des flottants.
De plus quand vous utilisez des flottants il y a souvent un moment où vous devez les convertir en entiers ou
l'inverse.

Avec les entiers, les additions et soustractions sont plus rapides que les multiplications et les divisions.
En fait, c'était vrai en 1962, maintenant c'est seulement partiellement vrai.
Les divisions sont toujours plus lentes (à moins que ça soit par une puissance de 2), mais les processeurs modernes
utilisent des algorithmes qui rendent les multiplications aussi rapides que des additions.
Quoi qu'il en soit, on va considérer que nous sommes toujours en 1962 parce que je veux expliquer l'algorithme
d'origine.
Avec les flottants, les multiplications sont plus lentes que les additions sur la plupart des processeurs modernes.

Soyez toujours prudents avec les fonctions que vous n'avez pas codé vous-même.
Est-ce que gfx.setPixel() est rapide ? Est-ce qu'elle prend toujours le même temps quand vous l'appelez 2 fois avec
les mêmes paramètres ?
Eh bien, comme je l'ai codée moi-même, je peux vous dire qu'elle n'est peut-être pas la fonction la plus efficace
pour cet algo particulier, mais elle devrait être assez rapide.
De toute façon on ne va pas la recoder parce que ce n'est pas le propos de cet article.
Le point important c'est que même si on optimise tout le reste dans cet algorithme, elle sera toujours appelée une
fois pour chaque pixel de la ligne. Donc sa durée totale d'exécution ne changera pas.

Alors maintenant, regardons les optimisations que Bresenham a faites les unes après les autres.

Supprimer la multiplication

Regardez la formule à l'intérieur de la boucle: 
		float y = slope * (x - x0) + y0;
Comment y change quand x va de x0 à x1 ?

Au début, quand x = x0, (x - x0) s'annule et y démarre à y0.

Ensuite, à chaque fois que x augmente de 1, y est augmenté de slope.

Donc on peut réécrire notre boucle: 
		float slope = (float)(y1 - y0) / (float)(x1 - x0);
		float y = y0;

		for (int x = x0; x <= x1; ++x)
		{
			gfx.setPixel(x, y, Color(255, 255, 255));
			y += slope;
		}
		gfx.render();
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
C'est une grande amélioration. Il y a beaucoup moins de calculs à l'intérieur de la boucle.
Il ne reste qu'une addition !

Maintenant, regardons le résultat:


Il y a toujours quelque chose qui ne va pas avec cette ligne.
En réalité c'est une erreur d'arrondi comme pour la première.
Mais vous pouvez vous demander où est l'erreur d'arrondi quand on n'utilise que des nombres en virgule flottante ?
C'est une des choses désagréables quand on utilise le C ou le C++: il fait des choses cachées.

Notre variable est un flottant. Le seul calcul qu'on fait dessus est une addition avec une autre variable flottante,
donc ça reste toujours un flottant.
Mais la fonction setPixel() prend des entiers en paramètres. Donc le compilateur C convertit implicitement notre y
flottant en entier.
Et quand il fait ça, il prend toujours l'entier inférieur.
Si notre y a une valeur de 5.7 par exemple, la conversion en entier donnera 5.
Mais on a vu au début que l'on devait prendre l'entier le plus proche pour avoir le pixel le plus proche de la
ligne idéale.

Alors on a vraiment besoin de corriger ce problème d'arrondi. Et la façon dont on va le faire est probablement la
partie la plus intéressante de l'algorithme de Bresenham.

Corriger l'arrondi

Comme on l'a dit, à chaque fois que x augmente de 1, y est augmenté de slope.
y doit être un flottant parce que la slope est un flottant.
Mais au final on affiche un pixel dont les coordonnées sont entières.
Alors séparons ces valeurs. D'un coté on aura un y entier qui représente la véritable coordonnée du pixel courant.
Et de l'autre on aura une valeur flottante qui contient la somme des slope (le reste).
Regardons en détail ce qui se passe dans cette image zoomée.


Ici on a la même pente que dans notre code: 10/30 = 1/3.
Remarquez que la ligne idéale qu'on veut tracer démarre au centre du premier pixel.

Pour le premier pixel la valeur de y est y0 et on dessine le pixel en y0.
Pour le 2ième pixel la valeur de y est y0 + 0.333... et on dessine encore le pixel en y0.
Pour le 3ième pixel la valeur de y est y0 + 0.666... et on dessine le pixel en (y0 + 1).
Remarquez qu'on a changé la coordonnée du pixel qu'on dessine quand la valeur flottante a dépassé 0.5.

Maintenant, nommons les variables pour clarifier les choses.
On a donc la relation:
		oldY = y + error
Maintenant récapitulons notre algorithme:
On commence par mettre y à y0 et error à 0
A chaque étape:
Mais pour respecter oldy = y + error, si on ajoute 1 à y, il faut aussi qu'on retire 1 à error.

Ca nous amène à ce code: 
		float slope = (float)(y1 - y0) / (float)(x1 - x0);
		int y = y0;
		float error = 0.0;

		for (int x = x0; x <= x1; ++x)
		{
			gfx.setPixel(x, y, Color(255, 255, 255));
			error += slope;

			if (error >= 0.5)
			{
				y++;
				error -= 1.0;
			}
		}
		gfx.render();

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Et voilà ce qu'on obtient:


Le fait de comparer error à 0.5 a corrigé notre problème d'arrondi.
Si vous prenez un screenshot de ce programme, que vous le chargez dans the Gimp et que vous dessinez une ligne
d'une autre couleur par dessus, vous verrez que vous obtenez exactement les mêmes pixels.
Parce the Gimp, comme la plupart des logiciels de dessin utilise l'algo de Bresenham lui aussi.

Supprimer les flottants

Continuons avec nos optimisations.
D'abord on va décaler l'erreur.
error démarre à 0.0 puis on la compare avec 0.5.
Ca serait exactement la meme chose si on la faisait démarrer à -0.5 puis qu'on la comparait à 0.0.
Comparer une valeur avec 0 est un peu plus rapide que de la comparer à une constante. Ca sera encore plus vrai
quand on l'aura convertie en entier.
Donc notre code devient:
		int y = y0;

		float slope = (float)(y1 - y0) / (float)(x1 - x0);
		float error = -0.5;

		for (int x = x0; x <= x1; ++x)
		{
			gfx.setPixel(x, y, Color(255, 255, 255));
			error += slope;

			if (error >= 0.0)
			{
				y++;
				error -= 1.0;
			}
		}
		gfx.render();

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Maintenant on va faire un petit travail de renommage pour que ça soit plus clair à l'étape suivante.
On va appeler "dx" la valeur (x1 - x0).
On va appeler "dy" la valeur (y1 - y0).
Et enfin on va appeler errorInc la valeur qu'on ajoute à error quand y change, c'est à dire -1.0.
		int y = y0;

		int dx = x1 - x0;
		int dy = y1 - y0;
		float slope = (float)dy / (float)dx;
		float error = -0.5;
		float errorInc = -1.0;

		for (int x = x0; x <= x1; ++x)
		{
			gfx.setPixel(x, y, Color(255, 255, 255));
			error += slope;

			if (error >= 0.0)
			{
				y++;
				error += errorInc;
			}
		}
		gfx.render();

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Maintenant on va convertir les valeurs flottantes en entiers.
Les seules variables flottantes qu'on a sont celles utilisées pour le calcul de error: slope, error et errorInc.
Pour éviter de perdre de la précision, l'idée est de multiplier ces variables par un nombre qui les rend toutes
entières.

slope vaut dy/dx donc on doit la multiplier au moins par dx pour la rendre entière.
error commence à -0.5 donc on doit la multiplier au moins par 2.
Finalement, si on multiplie toutes les variables par 2*dx elles deviendront toutes entières. 
		int y = y0;

		int dx = x1 - x0;
		int dy = y1 - y0;
		int slope = 2 * dy;
		int error = -dx;
		int errorInc = -2 * dx;

		for (int x = x0; x <= x1; ++x)
		{
			gfx.setPixel(x, y, Color(255, 255, 255));
			error += slope;

			if (error >= 0)
			{
				y++;
				error += errorInc;
			}
		}
		gfx.render();
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows

Est-ce que ça en valait la peine ?

Maintenant on a fait tout ce qu'on pouvait pour optimiser le code.
Arrivés à ce point, les articles qui expliquent l'algo de Bresenham essayent généralement de réduire le nombre de
variables ou d'inverser le signe d'error et d'errorInc.
Ce ne sont que des modifications esthétiques qui n'améliorent pas vraiment la vitesse du code.

Alors maintenant on va essayer de voir si toutes ces optimisations ont vraiment accéléré le code.
La lib SDL a un timer de haute précision qui est prévu justement pour ce genre de tâches.
J'ai ajouté 2 fonctions à ma lib dans la classe System: 
		void    StartPerfCounter();
		float   StopPerfCounter();
Vous avez juste à appeler StartPerfCounter() au début de la fonction que vous voulez chronométrer, et
StopPerfCounter() à la fin.
StopPerfCounter() renvoie un float qui est la nombre de millisecondes écoulées entre les deux appels.
J'ai utilisé ces fonctions pour chronométrer d'une part la première routine qu'on a écrite au début de l'article et
et d'autre part la dernière qu'on a écrite.
Pour avoir plus de précision, j'ai dessiné une série de lignes qui s'étendent sur toute le largeur de la fenêtre et
qui suivent les contraintes qu'on a données au départ.
Voici le code complet de ce test: 
		#include <stdio.h>
		#include <stdlib.h>
		#include <math.h>
		#include "main.h"
		#include "Graphics.h"
		#include "System.h"

		#define SCREEN_WIDTH    640
		#define SCREEN_HEIGHT   480

		void drawLineSlow(int x0, int y0, int x1, int y1, Color c)
		{
			float slope = (float)(y1 - y0) / (float)(x1 - x0);

			for (int x = x0; x <= x1; ++x)
			{
				float y = slope * (x - x0) + y0;
				gfx.setPixel(x, y, c);
			}
		}

		void drawLine(int x0, int y0, int x1, int y1, Color c)
		{
			int y = y0;

			int dx = x1 - x0;
			int dy = y1 - y0;
			int slope = 2 * dy;
			int error = -dx;
			int errorInc = -2 * dx;

			for (int x = x0; x <= x1; ++x)
			{
				gfx.setPixel(x, y, c);
				error += slope;

				if (error >= 0)
				{
					y++;
					error += errorInc;
				}
			}
		}

		int main(int argc, char* argv[])
		{
			// init the window
			gfx.init("Bresenham Line 7", SCREEN_WIDTH, SCREEN_HEIGHT);
			gfx.init2D();

			// test the slow function
			gfx.clearScreen(Color(0, 0, 0, SDL_ALPHA_OPAQUE));
			sys.StartPerfCounter();

			for (int i = 10; i < SCREEN_WIDTH; ++i)
				drawLineSlow(0, 0, i, 10, Color(i & 0xff, 255, 0));

			float time1 = sys.StopPerfCounter();
			printf("time1: %f ms\n", time1);

			// test the fast function
			gfx.clearScreen(Color(0, 0, 0, SDL_ALPHA_OPAQUE));
			sys.StartPerfCounter();

			for (int i = 10; i < SCREEN_WIDTH; ++i)
				drawLine(0, 0, i, 10, Color(i & 0xff, 255, 0));

			float time2 = sys.StopPerfCounter();
			printf("time2: %f ms\n", time2);

			gfx.render();

			while (sys.isQuitRequested() == false)
			{
				sys.processEvents();
			}

			gfx.quit();

			return EXIT_SUCCESS;
		}
		
		Télécharger le code source
		Télécharger l'exécutable pour Windows
En l'exécutant 3 fois sur mon ordinateur qui a un Intel Core i3-4330, j'obtiens une moyenne de 0.938 ms pour la
première routine et 0.514 ms pour celle qui est optimisée.
Donc, même avec un processeur moderne et avec les optimisations du compilateur on a écrit un code qui est environ
2 fois plus rapide.
C'est pas mal du tout.
Et rappelez-vous que ce temps inclut le temps d'exécution de setPixel() qu'on n'a pas optimisé.

Généraliser l'algorithme

Si on utilise notre code optimisé pour dessiner une ligne qui est presque verticale, comme 
		int x0 = 0;
		int y0 = 0;
		int x1 = 10;
		int y1 = 30;
On obtient une ligne à 45 degrés:


Avec la première fonction qu'on a écrite on obtient un autre résultat.


Il y a des pixels qui manquent.
C'est assez facile à comprendre car notre fonction ne dessine qu'un pixel pour chaque coordonnée x.
Pour des lignes qui sont plus verticales on devrait faire une boucle sur y au lieu de x.

Maintenant, rappelez-vous les contraintes qu'on a imposées au début: On peut diviser le plan de l'écran en 8 "octants".


Les contraintes ci-dessus nous limitent à l'octant 0.
A ce moment en général les autres articles se contentent d'écrire une routine séparée pour chaque octant et 2 de
plus pour les cas des lignes parfaitement horizontale ou verticale.
C'est dans l'esprit de l'algorithme pour avoir le code le plus rapide possible.
Mais je préfère un code plus compact même s'il est un peu plus lent. Ca sera beaucoup plus facile si on doit le
modifier dans le futur.

Observons les octants 0, 3, 4 et 7, c'est tous les cas où la ligne est plus horizontale que verticale.
Pour tous ces cas on va garder la boucle à travers les coordonnées x. Mais la direction peut changer: x1 peut être
inférieur à x0 ou y1 peut être inférieur à y0.

On va introduire 2 macros mathématiques.
ABS retourne la valeur absolue d'un nombre.
		#define ABS(_x) ((_x) >= 0 ? (_x) : -(_x))
SGN retourne le signe d'une valeur: soit -1, 1 ou 0 si le nombre est négatif, positif ou nul.
		#define SGN(_x) ((_x) < 0 ? -1 : \
		                 ((_x) > 0 ? 1 : 0))
Maintenant on va utiliser 2 variables d'incrément pour x et y pour gérer tous les cas. 
		void drawLine(int x0, int y0, int x1, int y1, Color c)
		{
			int y = y0;

			int dx = x1 - x0;
			int dy = y1 - y0;
			int incX = SGN(dx);
			int incY = SGN(dy);
			int slope = 2 * ABS(dy);
			int error = -ABS(dx);
			int errorInc = -2 * ABS(dx);

			for (int x = x0; x != x1 + incX; x += incX)
			{
				gfx.setPixel(x, y, c);
				error += slope;

				if (error >= 0)
				{
					y += incY;
					error += errorInc;
				}
			}
		}

		Télécharger le code source
		Télécharger l'exécutable pour Windows
Ensuite c'est facile d'écrire tous les cas. Les cas verticaux sont juste le symétrique des horizontaux. 
		#define ABS(_x) ((_x) >= 0 ? (_x) : -(_x))
		#define SGN(_x) ((_x) < 0 ? -1 : \
						 ((_x) > 0 ? 1 : 0))

		void drawLine(int x0, int y0, int x1, int y1, Color c)
		{
			int dx = x1 - x0;
			int dy = y1 - y0;
			int incX = SGN(dx);
			int incY = SGN(dy);
			dx = ABS(dx);
			dy = ABS(dy);

			if (dy == 0)
			{
				// horizontal line
				for (int x = x0; x != x1 + incX; x += incX)
					gfx.setPixel(x, y0, c);
			}
			else if (dx == 0)
			{
				// vertical line
				for (int y = y0; y != y1 + incY; y += incY)
					gfx.setPixel(x0, y, c);
			}
			else if (dx >= dy)
			{
				// more horizontal than vertical
				int slope = 2 * dy;
				int error = -dx;
				int errorInc = -2 * dx;
				int y = y0;

				for (int x = x0; x != x1 + incX; x += incX)
				{
					gfx.setPixel(x, y, c);
					error += slope;

					if (error >= 0)
					{
						y += incY;
						error += errorInc;
					}
				}
			}
			else
			{
				// more vertical than horizontal
				int slope = 2 * dx;
				int error = -dy;
				int errorInc = -2 * dy;
				int x = x0;

				for (int y = y0; y != y1 + incY; y += incY)
				{
					gfx.setPixel(x, y, c);
					error += slope;

					if (error >= 0)
					{
						x += incX;
						error += errorInc;
					}
				}
			}
		}

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

Test

Finalement, j'ai mis cette routine de ligne dans Graphics et j'ai écrit une sorte d'économiseur d'écran pour tester
des lignes dans toutes les orientations.


		#include <stdio.h>
		#include <stdlib.h>
		#include <math.h>
		#include <time.h>
		#include "main.h"
		#include "Graphics.h"
		#include "System.h"

		#define SCREEN_WIDTH    640
		#define SCREEN_HEIGHT   480

		int main(int argc, char* argv[])
		{
			// init the window
			gfx.init("Bresenham Line 10", SCREEN_WIDTH, SCREEN_HEIGHT);
			gfx.init2D();
			gfx.clearScreen(Color(0, 0, 0, SDL_ALPHA_OPAQUE));

			srand(time(NULL));

			int x0 = rand() % SCREEN_WIDTH;
			int y0 = rand() % SCREEN_HEIGHT;
			int spdX0 = (rand() % 7) - 3;
			int spdY0 = (rand() % 7) - 3;

			int x1 = rand() % SCREEN_WIDTH;
			int y1 = rand() % SCREEN_HEIGHT;
			int spdX1 = (rand() % 7) - 3;
			int spdY1 = (rand() % 7) - 3;

			int r = rand() % 256;
			int g = rand() % 256;
			int b = rand() % 256;
			int spdR = (rand() % 7) - 3;
			int spdG = (rand() % 7) - 3;
			int spdB = (rand() % 7) - 3;

			while (sys.isQuitRequested() == false)
			{
				gfx.line(x0, y0, x1, y1, Color(r, g, b));

				x0 += spdX0;
				if (x0 < 0 || x0 >= SCREEN_WIDTH)
				{
					spdX0 = -spdX0;
					x0 += spdX0;
				}

				y0 += spdY0;
				if (y0 < 0 || y0 >= SCREEN_HEIGHT)
				{
					spdY0 = -spdY0;
					y0 += spdY0;
				}

				x1 += spdX1;
				if (x1 < 0 || x1 >= SCREEN_WIDTH)
				{
					spdX1 = -spdX1;
					x1 += spdX1;
				}

				y1 += spdY1;
				if (y1 < 0 || y1 >= SCREEN_HEIGHT)
				{
					spdY1 = -spdY1;
					y1 += spdY1;
				}

				r += spdR;
				if (r < 0 || r > 255)
				{
					spdR = -spdR;
					r += spdR;
				}

				g += spdG;
				if (g < 0 || g > 255)
				{
					spdG = -spdG;
					g += spdG;
				}

				b += spdB;
				if (b < 0 || b > 255)
				{
					spdB = -spdB;
					b += spdB;
				}

				gfx.render();
				sys.wait(10);
				sys.processEvents();
			}

			gfx.quit();

			return EXIT_SUCCESS;
		}

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

Liens

Vidéo du dernier programme