diff --git a/README.md b/README.md new file mode 100644 index 0000000..bf4ac89 --- /dev/null +++ b/README.md @@ -0,0 +1,21 @@ +# Mandelbrot + +Esta es una colección de diferentes scripts para generar conjuntos de Mandelbrot. + +Implementada C, con aritmética flotante de 96 y 128 bits, paralelizada en OpenMP y con interfaz sencilla en SDL2. + +Los keyframes están generados en C mientras que las imágenes intermedias se interpolan con ImageMagik + +## Ejemplos: + +https://www.youtube.com/watch?v=j7vwzfUibwY + +https://www.youtube.com/watch?v=_8ECNKlCjbw + +https://www.youtube.com/watch?v=tcUKCw993Ms + +### [Buddahbrots](https://es.wikipedia.org/wiki/Buddhabrot) + +![](buddahbrot.png) + +![](buddahbrot2.png) diff --git a/buddahbrot.png b/buddahbrot.png new file mode 100644 index 0000000..c1a8ef2 Binary files /dev/null and b/buddahbrot.png differ diff --git a/buddahbrot2.png b/buddahbrot2.png new file mode 100644 index 0000000..0dbdd7b Binary files /dev/null and b/buddahbrot2.png differ diff --git a/generate.py b/generate.py new file mode 100644 index 0000000..4cc2411 --- /dev/null +++ b/generate.py @@ -0,0 +1,14 @@ +import os + +s = (1920, 1080) +f = 0.5 ** (1/69) + +j = 1 +for x in range(5): + w = 1920.0 + for i in range(69): + w /= f + width = int(w + 0.5) + print("convert img2-%05d.png -resize %dx1080^ -gravity Center -extent 1920x1080 img-%05d.png" % (x, width, 4160 + j)) + os.system("convert img2-%05d.png -resize %dx1080^ -gravity Center -extent 1920x1080 img-%05d.png" % (x, width, 4160 + j)) + j += 1 diff --git a/mandelbrot_buddah.c b/mandelbrot_buddah.c new file mode 100644 index 0000000..1339085 --- /dev/null +++ b/mandelbrot_buddah.c @@ -0,0 +1,118 @@ +#include +#include +#include +#include +#include + +#define VENTANA_ANCHO (1920*2) +#define VENTANA_ALTO (1080*2) +#define ANCHO VENTANA_ANCHO +#define ALTO VENTANA_ALTO +#define MAX_BAILOUT 50000 + +double xc = -0.5, yc = 0; + +int bailout(double x0, double y0, double zoom) { + x0 = x0 * zoom + xc; + y0 = y0 * zoom + yc; + + if( + (x0 > -1.2 && x0 <= -1.1 && y0 > -0.1 && y0 < 0.1) + || (x0 > -1.1 && x0 <= -0.9 && y0 > -0.2 && y0 < 0.2) + || (x0 > -0.9 && x0 <= -0.8 && y0 > -0.1 && y0 < 0.1) + || (x0 > -0.69 && x0 <= -0.61 && y0 > -0.2 && y0 < 0.2) + || (x0 > -0.61 && x0 <= -0.5 && y0 > -0.37 && y0 < 0.37) + || (x0 > -0.5 && x0 <= -0.39 && y0 > -0.48 && y0 < 0.48) + || (x0 > -0.39 && x0 <= 0.14 && y0 > -0.55 && y0 < 0.55) + || (x0 > 0.14 && x0 < 0.29 && y0 > -0.42 && y0 < -0.07) + || (x0 > 0.14 && x0 < 0.29 && y0 > 0.07 && y0 < 0.42) + ) return MAX_BAILOUT; + + double x = 0; + double y = 0; + int n = 0; + while(x*x + y*y <= 4 && n < MAX_BAILOUT) { + double xtemp = x*x - y*y + x0; + y = 2*x*y + y0; + x = xtemp; + n++; + } + + return n; +} + +unsigned int max(unsigned int m[ALTO][ANCHO]) { + unsigned int max = 0; + + for(int i = 0; i < ALTO; i++) + for(int j = 0; j < ANCHO; j++) + if(m[i][j] > max) + max = m[i][j]; + return max; +} + +void computar_intensidad(unsigned int m500[ALTO][ANCHO], unsigned int m5000[ALTO][ANCHO], unsigned int m50000[ALTO][ANCHO], double x0, double y0, double zoom) { + x0 = x0 * zoom + xc; + y0 = y0 * zoom + yc; + double x = 0; + double y = 0; + int n = 0; + while(x*x + y*y <= 4 && n < MAX_BAILOUT) { + double xtemp = x*x - y*y + x0; + y = 2*x*y + y0; + x = xtemp; + + int xx = (x - xc) / zoom + ANCHO / 2; + int yy = (y - yc) / zoom + ALTO / 2; + if(xx > 0 && xx < ANCHO && yy > 0 && yy < ALTO) { + if(n < 500) + m500[yy][xx]++; + if(n < 5000) + m5000[yy][xx]++; + m50000[yy][xx]++; + } + + n++; + } +} + + +int main() { + double zoom = 1e-3; + float step = 0.1; + + unsigned int (*m500)[ANCHO] = calloc(ANCHO * ALTO, sizeof(unsigned int)); + unsigned int (*m5000)[ANCHO] = calloc(ANCHO * ALTO, sizeof(unsigned int)); + unsigned int (*m50000)[ANCHO] = calloc(ANCHO * ALTO, sizeof(unsigned int)); + + for(double vy = -ALTO / 2; vy < ALTO / 2; vy += step) + for(double vx = - ANCHO / 2; vx < ANCHO / 2; vx += step) { + if(bailout(vx, vy, zoom) < MAX_BAILOUT) + computar_intensidad(m500, m5000, m50000, vx, vy, zoom); + } + + int max500 = max(m500); + int max5000 = max(m5000); + int max50000 = max(m50000); + +/* if(!max500) max500 = 1; + if(!max5000) max5000 = 1; + if(!max50000) max50000 = 1;*/ + int mm = max500; + if(max5000 > mm) + mm = max5000; + if(max50000 > mm) + mm = max50000; + + printf("P3\n%d %d\n%d\n", ALTO, ANCHO, mm); + for(int i = 0; i < ANCHO; i++) + for(int j = 0; j < ALTO; j++) + printf("%u %u %u\n", + m50000[j][i], //* 255 / max50000, + m5000[j][i], //* 255 / max5000, + m500[j][i] //* 255 / max500 + ); + + return 0; +} + diff --git a/mandelbrot_main.c b/mandelbrot_main.c index d0f5153..5591848 100644 --- a/mandelbrot_main.c +++ b/mandelbrot_main.c @@ -4,11 +4,11 @@ #include -#define VENTANA_ANCHO 300 -#define VENTANA_ALTO 200 +#define VENTANA_ANCHO 800 +#define VENTANA_ALTO 600 #define ANCHO VENTANA_ANCHO #define ALTO VENTANA_ALTO -#define MAX_ITER 10000 +#define MAX_ITER 500 #define JUEGO_FPS 10 static double time_in_ms(void) { @@ -22,8 +22,7 @@ static double time_in_ms(void) { //long double xc = -1.740119020442992652289958266376L; //long double yc = 0.028019209975637676263081582255L; -//long double xc = 0.264226442163892496529931097626331393L, yc = 0.002572261418806922467217290773078275L; // Video 2 -long double xc = -1.749900374662927212679117139337847675L, yc = 0.000000000002454985822810237255224056L; +long double xc = 0.264226442163892496529931097626331393L, yc = 0.002572261418806922467217290773078275; uint32_t computar_intensidad(long double x0, long double y0, long double zoom) { //x0 = x0 * zoom - 1.7400623825793398724575; // Derecha: Achicar @@ -46,7 +45,7 @@ uint32_t computar_intensidad(long double x0, long double y0, long double zoom) { return 0; //float h = (n + 1 - log(log2(sqrt(x*x+y*y)))) / MAX_ITER * 360; - float h = (int)(n * 360.0 / 100) % 360; + float h = (int)(n * 360.0 / 50) % 360; float xxx = (1 - fabs(fmodf(h / 60.0, 2) - 1)); uint8_t xx = (xxx < 1 ? xxx : 1) * 255; uint8_t r = 0, g = 0, b = 0; @@ -154,7 +153,7 @@ int width, height; for(int vx = - ANCHO / 2; vx < ANCHO / 2; vx++) { canvas[i++] = computar_intensidad(vx, vy, zoom); } -printf("%f %.5Le\nlong double xc = %.36LfL, yc = %.36LfL;\n", time_in_ms() - begin, zoom, xc, yc); +printf("%f %.5Le %.36Lf %.36Lf\n", time_in_ms() - begin, zoom, xc, yc); // Implementar imagen_a_textura() cuanto antes! diff --git a/mandelbrot_png.c b/mandelbrot_png.c index 5349624..74cca05 100644 --- a/mandelbrot_png.c +++ b/mandelbrot_png.c @@ -4,8 +4,6 @@ #include #include -#include - #define PREFIX "img-" #define VENTANA_ANCHO 1920 @@ -66,17 +64,14 @@ int main() { long double zoom = 0.179; double factor = 0.99; - size_t n = 0; - for(double z = zoom; z >= 1e-20; z *= factor) - n++; + int i = 0; - double *zooms = malloc(n * sizeof(double)); - zooms[0] = zoom; - for(size_t i = 1; i < n; i++) - zooms[i] = zooms[i - 1] * factor; + while(1) { + if(zoom < 1e-20) + break; + + zoom *= factor; - #pragma omp parallel for schedule(dynamic) - for(int i = 0; i < n; i++) { printf("%d\n", i); char aux[100]; @@ -87,13 +82,15 @@ int main() { for(int vy = ALTO / 2; vy > - ALTO / 2; vy--) for(int vx = - ANCHO / 2; vx < ANCHO / 2; vx++) { - uint32_t c = computar_intensidad(vx, vy, zooms[i]); + uint32_t c = computar_intensidad(vx, vy, zoom); fprintf(f, "%d %d %d\n", c >> 24, (c >> 16) & 0xFF, (c >> 8) & 0xFF); } fclose(f); sprintf(aux, "convert %s%05d.ppm %s%05d.png; rm %s%05d.ppm", PREFIX, i, PREFIX, i, PREFIX, i); system(aux); + + i++; } system("ffmpeg -y -framerate 24 -i " PREFIX "%05d.png -c:v libx264 -crf 0 output.mp4"); diff --git a/mandelbrot_png_openmp.c b/mandelbrot_png_openmp.c new file mode 100644 index 0000000..65cb1bd --- /dev/null +++ b/mandelbrot_png_openmp.c @@ -0,0 +1,114 @@ +#include +#include +#include +#include +#include + +#include + +#include + +#define PREFIX "img2-" + +#define VENTANA_ANCHO (1920*2) +#define VENTANA_ALTO (1080*2) +#define ANCHO VENTANA_ANCHO +#define ALTO VENTANA_ALTO +#define MAX_ITER 10000 //500 en video2 + +//__float128 xc = -1.740119020442992652289958266376L; +//__float128 yc = 0.028019209975637676263081582255L; +//__float128 xc = 0.264226442163892496529931097626331393L, yc = 0.002572261418806922467217290773078275; +__float128 xc = -1.749900374662927212679117139337847675Q, yc = 0.000000000002454985822810237255224056Q; + +uint32_t computar_intensidad(__float128 x0, __float128 y0, __float128 zoom) { + //x0 = x0 * zoom - 1.7400623825793398724575; // Derecha: Achicar + //y0 = y0 * zoom + 0.0281753397792111; // Abajo: Achicar + x0 = x0 * zoom + xc; // Derecha: Achicar + y0 = y0 * zoom + yc; // Abajo: Achicar + //x0 = x0 * zoom - 1.78; // Derecha: Achicar + //y0 = y0 * zoom + 0.0; // Abajo: Achicar + __float128 x = 0; + __float128 y = 0; + int n = 0; + while(x*x + y*y <= 4 && n < MAX_ITER) { + __float128 xtemp = x*x - y*y + x0; + y = 2*x*y + y0; + x = xtemp; + n++; + } + + if(n == MAX_ITER) + return 0; + + //float h = (n + 1 - log(log2(sqrt(x*x+y*y)))) / MAX_ITER * 360; + //float h = n * 360.0 / MAX_ITER; + float h = (int)(n * 360.0 / 100) % 360; + float xxx = (1 - fabs(fmodf(h / 60.0, 2) - 1)); + uint8_t xx = (xxx < 1 ? xxx : 1) * 255; + uint8_t r = 0, g = 0, b = 0; + + if(h < 60) + r = 255, g = xx; + else if(h < 120) + g = 255, r = xx; + else if(h < 180) + g = 255, b = xx; + else if(h < 240) + b = 255, g = xx; + else if(h < 300) + b = 255, r = xx; + else + r = 255, b = xx; + + return (r << 24) | (g << 16) | (b << 8); +} + + +int main() { + __float128 zoom = 0.179; + double factor = 0.99; + + size_t n = 0; + for(double z = zoom; z >= 1e-21; z *= factor) + n++; + + double *zooms = malloc(n * sizeof(double)); + zooms[0] = zoom; + for(size_t i = 1; i < n; i++) + zooms[i] = zooms[i - 1] * factor; + + zooms[0] = zooms[4160] * 0.5; + printf("%.16e\n", zooms[4160]); + n = 6; + for(size_t i = 1; i < n; i++) + zooms[i] = zooms[i - 1] * 0.5; + + #pragma omp parallel for schedule(dynamic) + for(int i = 0; i < n; i++) { + printf("%d %.16e\n", i, zooms[i]); + + char aux[100]; + sprintf(aux, "%s%05d.ppm", PREFIX, i); + FILE *f = fopen(aux, "wt"); + fprintf(f, "P3\n%d %d\n255\n", ANCHO, ALTO); + + + for(int vy = ALTO / 2; vy > - ALTO / 2; vy--) + for(int vx = - ANCHO / 2; vx < ANCHO / 2; vx++) { + uint32_t c = computar_intensidad(vx, vy, zooms[i]); + fprintf(f, "%d %d %d\n", c >> 24, (c >> 16) & 0xFF, (c >> 8) & 0xFF); + } + fclose(f); + + sprintf(aux, "convert %s%05d.ppm %s%05d.png; rm %s%05d.ppm", PREFIX, i, PREFIX, i, PREFIX, i); + system(aux); + } + + free(zooms); + +// system("ffmpeg -y -framerate 24 -i " PREFIX "%05d.png -c:v libx264 -crf 0 output.mp4"); + + return 0; +} +