HALLEY MAP

/*

Calcul de la carte de Halley
 
L'ecran correspond a une portion du plan complexe
A chaque pixel (i, j) correspond donc un complexe z0

On considere la fonction H(z)= z- F(z)/ ( F'(z)- F"(z)* F(z) / ( 2* F'(z) ) ) 
 avec F une fonction du plan complexe admettant plusieurs zeros
Cette formule provient du developpement limite a l'ordre 2 de F
 
On itere H avec z0 comme valeur initiale : z1= H(z0); z2= H(z1); ... 
et on recupere le n minimum tel que module_carre( z(n+1)- zn )< EPS (pour epsilon) fixe
On considere que la suite diverge si n trop grand (> STOP)

On obtient ainsi une correspondance pixel <-> entier qui va fournir une image en niveaux de gris
Si pour un pixel donné, la suite diverge, on le remplit par un pixel noir. 

 
 
Precisions
 
 ici, j'ai pris comme exemples : F(z)= sin(z) ; F(z)= z* ( z^6- 1)
 Il est possible de faire des zooms successifs avec le clic gauche de la souris; 
 un clic droit fait revenir a la fenetre d'origine : [-2.5, 2.5] * [-2.5, 2.5]
 STOP== 256 de facon a avoir des niveaux de gris sur 8 bits
 
 
Problemes et ameliorations
 
 Plus on zoome, plus le programme est lent, je ne sais pas pourquoi; 
 il y a de toute facon un probleme de precision des double lorsque le zoom est trop important.
 
 Je multiplie par 10 la valeur des indices n, avant d'en faire un niveau de gris, sinon l'image est trop sombre
 Mais je tronque a 255, d'ou une perte d'info; il faudrait un reetalage de l'histogramme
 
 Une interface de choix de la fonction F serait mieux que changer une ligne de code !
*/


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <complex.h>

#include "SDL.h"

// voir commentaire ci-dessus
#define EPS 0.0001
#define STOP 256

// fenetre initiale
#define MIN_X_INIT -2.5
#define MAX_X_INIT 2.5
#define MIN_Y_INIT -2.5
#define MAX_Y_INIT 2.5

// nombre pixels
#define NPIX_X 800
#define NPIX_Y 800


// vont remplacer les constantes de fenetre initiale lors des zooms
double min_x, max_x, min_y, max_y;



inline int min(int a, int b)
{
	return a < b ? a : b;
}

inline int max(int a, int b)
{
	return a > b ? a : b;
}


// place la couleur pixel a l'emplacement (x, y) de la surface "surface" (voir doc SDL)
void putpixel(SDL_Surface *surface, int x, int y, Uint32 pixel)
{
    int bpp = surface->format->BytesPerPixel;
    Uint8 *p = (Uint8 *)surface->pixels + y * surface->pitch + x * bpp;
	
    switch(bpp) {
		case 1:
			*p = pixel;
			break;
			
		case 2:
			*(Uint16 *)p = pixel;
			break;
			
		case 3:
			if(SDL_BYTEORDER == SDL_BIG_ENDIAN) {
				p[0] = (pixel >> 16) & 0xff;
				p[1] = (pixel >> 8) & 0xff;
				p[2] = pixel & 0xff;
			} else {
				p[0] = pixel & 0xff;
				p[1] = (pixel >> 8) & 0xff;
				p[2] = (pixel >> 16) & 0xff;
			}
			break;
			
		case 4:
			*(Uint32 *)p = pixel;
			break;
    }
}


// Trace une ligne entre (x1, y1) et (x2, y2); pas optimise
void draw_line(SDL_Surface *surface, int x1, int y1, int x2, int y2, Uint32 color)
{
	int i, steep, tmp, dx, dy, error, ystep;
	
	if (abs(y2-y1) > abs(x2-x1)) steep= 1; else steep= 0;
	if (steep)
	{
		tmp= x1; x1= y1; y1= tmp;
		tmp= x2; x2= y2; y2= tmp;
	}
	if (x1> x2)
	{
		tmp= x1; x1= x2; x2= tmp;
		tmp= y1; y1= y2; y2= tmp;
	}
	dx= x2- x1;
	dy= abs(y2- y1);
	error= dx/ 2;
	tmp= y1;
	if (y1< y2) ystep= 1; else ystep= -1;
	for (i= x1; i<= x2; i++)
	{
		if (steep) putpixel(surface, tmp, i, color); else putpixel(surface, i, tmp, color);
		error-= dy;
		if (error< 0)
		{
			tmp+= ystep;
			error+= dx;
		}
	}
}


// Fonctions tests
complex f_test_1(complex z)
{
	return ( cpow(z, 7.)- z )/ ( 7.* cpow(z, 6.)- 1.- 42.* cpow(z, 5.)* (cpow(z, 7.)- z) / (14.* cpow(z, 6.)- 2.) );
}


complex f_test_2(complex z)
{
	return csin(2.* z)/( 1.+ ccos(z)* ccos(z) );
}


// renvoie l'indice du plus petit n tel que diff< EPS (voir commentaire de debut de page)
int cvgence(complex z)
{
	int n= 0;
	
	// p_f pointe sur une des fonctions tests correspondant a F
	complex (*p_f)(complex);
	// Il modifier cette ligne pour changer F
	// p_f= &f_test_1;
	p_f= &f_test_2;
	
	complex z_;
	
	while (1)
	{
		n++;
		z_= z- (*p_f)(z);
		if (cabs(z_- z)* cabs(z_- z)< EPS)
			break;
		if (n>= STOP)
			break;
		z= z_;
	}
	if (n>= STOP)
		// renvoie -1 si divergence
		return -1;
	else
		return n;
}


// Place dans le buffer buf les valeurs entre 0 et 255 correspondant aux niveaux de gris
// La carte de Halley est evaluee sur [min_x_, max_x_]* [min_y_, max_y_],
// sachant que NPIX_X et NPIX_Y, donc la taille de l'ecran, ne sont pas modifiables
void compute_halley(SDL_Surface *sdl_surf, Uint32 * buf, double min_x_, double max_x_, double min_y_, double max_y_)
{
	Uint32 color;
	double pas_x, pas_y;
	int i, j, n;
	complex z;
	
	// (pas_x, pas_y) est la taille d'un pixel dans le plan complexe
	pas_x= (max_x_- min_x_)/ (double)(NPIX_X);
	pas_y= (max_y_- min_y_)/ (double)(NPIX_Y);
	
	for (j=0; j<NPIX_Y; j++)
	{
		for (i=0; i<NPIX_X; i++)
		{
			z= ((double)(min_x_)+ (double)(i)* pas_x)+ I* ((double)(min_y_)+ (double)(j)* pas_y);
			n= cvgence(z);
			if (n< 0) /* On associe la valeur 0 (noir) a un complexe qui fait diverger la suite */
			{
				color= SDL_MapRGB(sdl_surf->format, 0, 0, 0);
			}
			else
			{
				// pour eclaircir l'image, multiplie par 10, et tronque a 255 (donc saturation)
				n*= 10;
				if (n> 256)
					n= 255;
				color= SDL_MapRGB(sdl_surf->format, n, n, n);
			}
			// remplir le buffer
			buf[i+ j* NPIX_X]= color;
		}
	}
	
}


// dessine le contenu du buffer rempli par compute_halley
void draw_halley(SDL_Surface *sdl_surf, Uint32 * buf)
{
	int i, j;
	for (j=0; j<NPIX_Y; j++)
	{
		for (i=0; i<NPIX_X; i++)
		{
			putpixel(sdl_surf, i, j, buf[i+ j* NPIX_X]);
		}
	}
}


// zoom sur le plus petit carré englobant le rectangle [x0_, x1_] * [y0_, y1_]
// de facon a eviter les distortions (applatissement et etirement)
void zoom(SDL_Surface *sdl_surf, Uint32 * buf, int x0_, int y0_, int x1_, int y1_)
{
	double pas_x, pas_y;
	double min_x_new, max_x_new, min_y_new, max_y_new;
	
	// taille pixel dans le plan complexe
	pas_x= (max_x- min_x)/ (double)(NPIX_X);
	pas_y= (max_y- min_y)/ (double)(NPIX_Y);
	
	// coordonnees dans le plan complexe des extremites du rectangle de selection
	min_x_new= min_x+ pas_x* (double)(min(x0_, x1_));
	max_x_new= min_x+ pas_x* (double)(max(x0_, x1_));
	min_y_new= min_y+ pas_y* (double)(min(y0_, y1_));
	max_y_new= min_y+ pas_y* (double)(max(y0_, y1_));
	
	// passage au format carré
	if (max_x_new- min_x_new> max_y_new- min_y_new)
	{
		min_y_new= (max_y_new+ min_y_new)/ 2. - (max_x_new- min_x_new)/ 2.;
		max_y_new= (max_y_new+ min_y_new)/ 2. + (max_x_new- min_x_new)/ 2.;
	}
	else
	{
		min_x_new= (max_x_new+ min_x_new)/ 2. - (max_y_new- min_y_new)/ 2.;
		max_x_new= (max_x_new+ min_x_new)/ 2. + (max_y_new- min_y_new)/ 2.;
	}
	
	// affectation des nouvelles valeurs
	min_x= min_x_new;
	max_x= max_x_new;
	min_y= min_y_new;
	max_y= max_y_new;
	
	compute_halley(sdl_surf, buf, min_x, max_x, min_y, max_y);
	draw_halley(sdl_surf, buf);
	
}



int main(int argc, char *argv[])
{
	Uint32 initflags = SDL_INIT_VIDEO;
	SDL_Surface *screen;
	Uint8  video_bpp = 0;
	// SDL_HWSURFACE pour la pris en charge du rendu par la carte graphique
	// SDL_DOUBLEBUF : plus rapide (voir doc SDL)
	Uint32 videoflags = SDL_HWSURFACE | SDL_DOUBLEBUF;

	// pour fermeture programme
	int done;
    
	// gestion des evenements SDL (souris, clavier)
	SDL_Event event;
	
	// booleen : vrai si bouton souris gauche enfonce
	int left_button_down= 0;
	
	// position du curseur lors du clic gauche (sert a selection zoom)
	int x0, y0;
	
	// couleur du rectangle de selection
	Uint32 color_selec;
	
	// buffer qui contiendra les valeurs pixels a afficher
	Uint32 *buf_halley;
	
	/* Initialiser SDL */
	if ( SDL_Init(initflags) < 0 ) {
		fprintf(stderr, "Couldn't initialize SDL: %s\n", SDL_GetError());
		exit(1);
	}

	/* Set video mode */
	screen=SDL_SetVideoMode(NPIX_X, NPIX_Y, video_bpp, videoflags);
        if (screen == NULL) {
		fprintf(stderr, "Couldn't set %ix%ix%i video mode: %s\n", video_bpp, NPIX_X, NPIX_Y, SDL_GetError());
		SDL_Quit();
		exit(2);
	}
	
	/* Initialisation des variables */
	done = 0;
	color_selec= SDL_MapRGB(screen->format, 100, 100, 0);
	
	buf_halley= malloc(sizeof(Uint32)* NPIX_X* NPIX_Y);

	min_x= MIN_X_INIT;
	max_x= MAX_X_INIT;
	min_y= MIN_Y_INIT;
	max_y= MAX_Y_INIT;
	
	// 1er calcul et affichage
	compute_halley(screen, buf_halley, MIN_X_INIT, MAX_X_INIT, MIN_Y_INIT, MAX_Y_INIT);
	draw_halley(screen, buf_halley);
	SDL_Flip(screen);
	
	while ( !done ) {

		while ( SDL_PollEvent(&event) ) {
			switch (event.type) {
				
				// mouvement souris : si bouton gauche enfonce, actualiser cadre de selection et rafraichir ecran
				case SDL_MOUSEMOTION:
					if (left_button_down)
					{
						draw_halley(screen, buf_halley);

						draw_line(screen, x0, y0, x0, event.motion.y, color_selec);
						draw_line(screen, event.motion.x, y0, event.motion.x, event.motion.y, color_selec);
						draw_line(screen, x0, y0, event.motion.x, y0, color_selec);
						draw_line(screen, x0, event.motion.y, event.motion.x, event.motion.y, color_selec);
						SDL_Flip(screen);
					}
					break;
				case SDL_MOUSEBUTTONDOWN:
					switch(event.button.button) {
						case SDL_BUTTON_LEFT:
							// si clic gauche, recuperer position curseur
							left_button_down= 1;
							x0= event.button.x;
							y0= event.button.y;
							break;
						
						case SDL_BUTTON_RIGHT:
							// si clic droit, on revient a l'emprise initiale
							min_x= MIN_X_INIT;
							max_x= MAX_X_INIT;
							min_y= MIN_Y_INIT;
							max_y= MAX_Y_INIT;
							
							compute_halley(screen, buf_halley, min_x, max_x, min_y, max_y);
							draw_halley(screen, buf_halley);
							SDL_Flip(screen);
							
							break;
						}
					break;
				case SDL_MOUSEBUTTONUP:
					switch(event.button.button) {
						// relache bouton gauche => zoom sur le cadre de selection
						case SDL_BUTTON_LEFT:
							left_button_down= 0;
							
							zoom(screen, buf_halley, x0, y0, event.button.x, event.button.y);
							SDL_Flip(screen);
							
							break;
						
						case SDL_BUTTON_RIGHT:
							break;
						}
					break;
				case SDL_KEYDOWN:
					// touche clavier => sortie programme
				case SDL_QUIT:
					done = 1;
					break;
				default:
					break;
			}
		}
	}
	
	/* Nettoyage */
	free(buf_halley);
	SDL_Quit();
	return(0);
}