Calcul de la carte de Halley et visualisation avec la librairie SDL
Génère une image en niveaux de gris de la carte de Halley (j'ai lu cet algo dans un livre de Pickover)
C'est une fractale du type Julia; a chaque complexe est associé une valeur dépendant de la vitesse de convergence d'une fonction a fournir (voir code)
Il est possible de zoomer en sélectionnant une partie de l'image
Compilé avec SDL
Source / Exemple :
/*
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:
break;
case 2:
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:
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);
}
Conclusion :
Largement améliorable j'imagine... Toute remarque est bienvenue
Vous n'êtes pas encore membre ?
inscrivez-vous, c'est gratuit et ça prend moins d'une minute !
Les membres obtiennent plus de réponses que les utilisateurs anonymes.
Le fait d'être membre vous permet d'avoir un suivi détaillé de vos demandes et codes sources.
Le fait d'être membre vous permet d'avoir des options supplémentaires.