385 lines
7.8 KiB
C
385 lines
7.8 KiB
C
/*
|
|
fractales.c
|
|
===========
|
|
|
|
*....................................................*
|
|
. warning, this is a 'work in progress' .
|
|
. .
|
|
. don't trust prototypes for building .
|
|
. serious applications... .
|
|
*....................................................*
|
|
|
|
*/
|
|
|
|
#include <stdio.h>
|
|
#include <stdlib.h>
|
|
#include <math.h>
|
|
#include <time.h>
|
|
|
|
#include <tthimage.h>
|
|
|
|
#include "fractales.h"
|
|
|
|
/*::------------------------------------------------------------------::*/
|
|
void Image_fractales_print_version(int flag)
|
|
{
|
|
printf("*** Fractales (%s) v 0.0.12\n", __FILE__);
|
|
if (flag)
|
|
{
|
|
printf("\tcompiled %s , %s\n", __DATE__, __TIME__);
|
|
}
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
/*
|
|
* fonction de comparaison de deux variables 'double'
|
|
*/
|
|
int compare_double(double a, double b)
|
|
{
|
|
|
|
return 0;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
/*
|
|
* cette fonction a été repompé (sans honte :) du livre de
|
|
* Roger T. Stevens, "Fractal Programming in C" éditée par
|
|
* M&T Books.
|
|
* Ceci dit, elle est un peu trop monolithique, puisque l'équation
|
|
* à résoudre est codée en dur.
|
|
*/
|
|
int Newton_0(Image_Desc *im, int maxiter)
|
|
{
|
|
int col, row, i;
|
|
double dcol, drow;
|
|
double Xmax, Xmin, Ymax, Ymin;
|
|
double deltaX, deltaY, denom;
|
|
double X, Y, Xsquare, Ysquare, Xold, Yold;
|
|
|
|
Xmax = Ymax = 3.5; Xmin = Ymin = -3.5;
|
|
deltaX = (Xmax-Xmin)/(double)im->width;
|
|
deltaY = (Ymax-Ymin)/(double)im->height;
|
|
|
|
for (col=0; col<im->width; col++)
|
|
{
|
|
dcol = (double)col;
|
|
for (row=0; row<im->height; row++)
|
|
{
|
|
drow = (double)row;
|
|
X = Xmin + dcol * deltaX;
|
|
Y = Ymin + drow * deltaY;
|
|
Xsquare = 0.0;
|
|
Ysquare = 0.0;
|
|
Xold = Yold = 42.0;
|
|
for (i=0; i<maxiter; i++)
|
|
{
|
|
Xsquare = X * X;
|
|
Ysquare = Y * Y;
|
|
denom = 3*((Xsquare-Ysquare)*(Xsquare-Ysquare) +
|
|
4*Xsquare*Ysquare);
|
|
if (denom < 0.0000001)
|
|
denom = 0.0000001;
|
|
|
|
X = .6666667*X + (Xsquare - Ysquare)/denom;
|
|
Y = .6666667*Y + 2*X*Y/denom;
|
|
if ((Xold==X) && (Yold==Y))
|
|
break;
|
|
} /* finbo sur maxiter */
|
|
|
|
if (X > 0)
|
|
{
|
|
Image_plotRGB(im, col, row, 0, 0, 0);
|
|
}
|
|
else
|
|
{
|
|
if (Y > 0)
|
|
Image_plotRGB(im, col, row, 255, 0, 0);
|
|
else
|
|
Image_plotRGB(im, col, row, 0, 0, 255);
|
|
}
|
|
}
|
|
if (0 == col%5)
|
|
fprintf(stderr, "newton 0: col %5d / %5d\r", col, im->width);
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
/*
|
|
* C'est pas encore fini ce truc ?
|
|
*/
|
|
int GingerBreadMan(Image_Desc *im, int maxiter)
|
|
{
|
|
double x1, y1, x2, y2;
|
|
int iter;
|
|
|
|
x1 = 10.0;
|
|
y1 = 0.333;
|
|
|
|
for (iter=0; iter<maxiter; iter++)
|
|
{
|
|
x2 = 1.0 - y1 + abs(x1);
|
|
y2 = x1;
|
|
|
|
printf("%9g %9g\n", x2, y2);
|
|
|
|
x1 = x2;
|
|
y1 = y2;
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
/*
|
|
En fait, il faudrait travailler sur deux images...
|
|
*/
|
|
int Image_Random_Walker_0(Image_Desc *im, Image_Desc *dst,
|
|
int nbpass, int nbsteps)
|
|
{
|
|
int pass, step, foo;
|
|
int x, y, ra, ga, ba, rb, gb, bb;
|
|
|
|
if ( (foo=Image_compare_desc(im, dst)) )
|
|
{
|
|
fprintf(stderr, "Random_Walker_0: images are differents %d\n", foo);
|
|
return foo;
|
|
}
|
|
|
|
for (pass=0; pass<nbpass; pass++)
|
|
{
|
|
x = rand() % im->width;
|
|
y = rand() % im->height;
|
|
|
|
ra = Image_R_pixel(im, x, y);
|
|
ga = Image_G_pixel(im, x, y);
|
|
ba = Image_B_pixel(im, x, y);
|
|
|
|
/* printf("\tP%d\t\tX%d\tY%d\t%d,%d,%d\n", pass, x, y, ra, ga, ba); */
|
|
|
|
for (step=0; step<nbsteps; step++)
|
|
{
|
|
x += (rand()%3)-1;
|
|
if (x<0) x += im->width;
|
|
if (x>=im->width) x -= im->width;
|
|
y += (rand()%3)-1;
|
|
if (y<0) y += im->height;
|
|
if (y>=im->height) y -= im->height;
|
|
|
|
rb = Image_R_pixel(im, x, y);
|
|
gb = Image_G_pixel(im, x, y);
|
|
bb = Image_B_pixel(im, x, y);
|
|
|
|
Image_plotRGB(dst, x, y,
|
|
(ra+rb+rb)/3, (ga+gb+gb)/3, (ba+bb+bb)/3);
|
|
}
|
|
}
|
|
|
|
return 0;
|
|
}
|
|
|
|
/*::------------------------------------------------------------------::*/
|
|
|
|
#define XMIN (-1.5)
|
|
#define XMAX (1.5)
|
|
#define YMIN (-1.5)
|
|
#define YMAX (1.5)
|
|
|
|
int Image_Julia_0_0(Image_Desc *im, RGB_map *pmap, int fuzz,
|
|
double cx, double cy)
|
|
{
|
|
int col, row;
|
|
double deltaX, deltaY, dx, dy, x2, y2;
|
|
long debut, fin;
|
|
int iter, r, g, b, maxiter;
|
|
|
|
debut = time(NULL);
|
|
|
|
fprintf(stderr, "Julia_0_0: %dx%d\n", im->width, im->height);
|
|
|
|
deltaX = (XMAX - XMIN) / im->width;
|
|
deltaY = (YMAX - YMIN) / im->height;
|
|
|
|
maxiter = 0;
|
|
|
|
/*
|
|
* si pas de palette définie, en créer une ?
|
|
*/
|
|
for (row=0; row<im->height; row++)
|
|
{
|
|
dy = YMIN + (double)row * deltaY;
|
|
|
|
for (col=0; col<im->width; col++)
|
|
{
|
|
#if DEBUG_LEVEL
|
|
printf("%5d %5d ", row, col);
|
|
#endif
|
|
|
|
dx = XMIN + (double)col * deltaX;
|
|
|
|
iter = 0;
|
|
x2 = y2 = 0.0;
|
|
|
|
while ( (iter < fuzz) &&
|
|
((x2 + y2) < 4.0) )
|
|
{
|
|
x2 = dx * dx;
|
|
y2 = dy * dy;
|
|
dx = 2 * dx * dy + cx;
|
|
dy = x2 - y2 + cy;
|
|
iter++;
|
|
}
|
|
|
|
if (iter > maxiter)
|
|
maxiter = iter;
|
|
|
|
if (iter >= fuzz)
|
|
{
|
|
r = g = b = 255;
|
|
}
|
|
else
|
|
{
|
|
r = pmap->red[iter%pmap->nbre];
|
|
g = pmap->green[iter%pmap->nbre];
|
|
b = pmap->blue[iter%pmap->nbre];
|
|
}
|
|
|
|
Image_plotRGB(im, col, row, r, g, b);
|
|
|
|
#if DEBUG_LEVEL
|
|
printf(" %3d \r", iter);
|
|
#endif
|
|
}
|
|
}
|
|
|
|
fin = time(NULL);
|
|
|
|
printf("duree du calcul: %ld secondes, maxiter: %d\n", fin-debut, maxiter);
|
|
|
|
return 0;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
|
|
#define XA -2.0
|
|
#define XB 2.0
|
|
#define YA -2.0
|
|
#define YB 1.3
|
|
|
|
int Image_Mandelbrot_0(Image_Desc *image, int maxiter)
|
|
{
|
|
int Y, X;
|
|
int r, g, b;
|
|
double zr, zi, fx, fy, real, imag, fiter, dist2;
|
|
long debut, fin, temps; /* pour chronometrer */
|
|
int iter, out, foo;
|
|
double tmpr, tmpi;
|
|
|
|
debut = time(NULL);
|
|
|
|
out = 0;
|
|
|
|
for (Y=0; Y<image->height; Y++)
|
|
{
|
|
fprintf(stderr, "ligne %4d sur %d\r", Y, image->height);
|
|
fy = (double)Y;
|
|
real = ( fy / (double)image->height * (YB-YA) ) + YA;
|
|
|
|
for (X=0; X<image->width; X++)
|
|
{
|
|
fx = (double)X;
|
|
imag = ( fx / (double)image->width * (XB-XA) ) + XA;
|
|
|
|
iter = 0;
|
|
zr = zi = 0.0;
|
|
|
|
do
|
|
{
|
|
iter++;
|
|
|
|
dist2 = zr*zr + zi*zi;
|
|
if (dist2 > 4.0)
|
|
{
|
|
out++;
|
|
break;
|
|
}
|
|
|
|
/*
|
|
* calculer le mandelbrot
|
|
* ----------------------
|
|
* la formule est:
|
|
* z0 = 0 + 0i
|
|
* z1 = z0^2 + c;
|
|
* ...
|
|
*/
|
|
tmpr = zr*zr - zi*zi;
|
|
tmpi = 2.0 * zr * zi;
|
|
|
|
zr = tmpr + real;
|
|
zi = tmpi + imag;
|
|
|
|
} while (iter < maxiter);
|
|
|
|
fiter = (double)iter / (double)maxiter;
|
|
|
|
Image_color_x(zr, zi, &r, &g, &b);
|
|
Image_plotRGB(image, X, Y, r, g, b);
|
|
|
|
}
|
|
|
|
fflush(stderr);
|
|
}
|
|
|
|
fin = time(NULL);
|
|
|
|
temps = fin - debut;
|
|
if (temps < 500)
|
|
{
|
|
foo = 's';
|
|
}
|
|
else
|
|
{
|
|
foo = 'm';
|
|
temps /= 60;
|
|
}
|
|
|
|
printf("\nMandelbrot_0: temps d'execution = %ld %c\n", temps, foo);
|
|
printf("Mandelbrot_0: nombre de points dehors = %d\n", out);
|
|
|
|
return 0;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
int Lorenz_Orbits(Image_Desc *img, int iters,
|
|
double a, double b, double c, double dt)
|
|
{
|
|
int pass;
|
|
double dx, dy, dz, dxp, dyp, dzp;
|
|
FILE *fp;
|
|
|
|
if (NULL == (fp=fopen("lorenz.dat", "w"))) {
|
|
perror("lorenz datafile");
|
|
exit(1);
|
|
}
|
|
|
|
fprintf(stderr, "%s ( %p %d / %f %f %f %f )\n", __func__, img, iters,
|
|
a, b, c, dt);
|
|
|
|
dx = dy = dz = 0.0001;
|
|
|
|
for (pass=0; pass<iters; pass++) {
|
|
|
|
fprintf(fp, "%6d %8f %8f %8f\n", pass, dx, dy, dz);
|
|
|
|
dxp = dx - (a*dx + a*dy) * dt;
|
|
dyp = dy + (b*dx - dy - dz*dz) * dt;
|
|
dzp = dz - (c*dz + dx*dy) * dt;
|
|
|
|
dx = dxp, dy = dyp, dz = dzp;
|
|
}
|
|
|
|
fclose(fp);
|
|
|
|
return 42;
|
|
}
|
|
/*::------------------------------------------------------------------::*/
|
|
/*
|
|
* et un PopCorn ?
|
|
*/
|
|
/*::------------------------------------------------------------------::*/
|