- Messaggi
- 2,332
- Reazioni
- 1,928
- Punteggio
- 134
Frattali: insieme di Mandelbrot
Questa volta voglio parlare di qualcosa di insolito rispetto al mio target. Anni fa scrissi un piccolissimo generatore di frattali, che pubblicai più avanti su Github, appartenendi all'insieme di Mandelbrot. Nulla di particolarmente avanzato e senza l'utilizzo di tecniche particolari di colorazione, sia chiaro; circa 1 mese fa ho aggiunto il supporto al multithreading.
Ho deciso di scrivere un brevissimo articolo solo per incuriosire chi magari non programma da molto tempo e non conosce i frattali. Spero anche di suscitare l'interese dei i più navigati, ovviamente!
Si trattava di una partecipazione ad un topic su un forum quella che mi fece scrivere queste due righe, quindi decisi come altri di avere poche dipendenze. Le immagini vengono generate in PPM; quelle che vedete nell'articolo sono quindi state tutte convertite.
L'articolo avrà il seguente svolgimento:
Immagini PPM: cosa sono?
Portable pixmap format, conosciuto anche come PPM, è un formato di immagine definito nel progetto Netpbm: allo stesso progetto appartengono PGM e PBM rispettivamente immagini in scala di grigi e in bianco e nero. PPM invece permette di rappresentare immagini a colori (0-255) in formato RGB. Sono molto semplici da gestire e possono essere in ASCII o binarie.
quello riportato qui sopra è l'header di un'immagine PPM rappresentata in binario (P6) larga w pixel e alta h pixel; ogni pixel dell'immagine ha una profondità di 24bit (8 * 3, dove 3 sono i canali di colore, RGB).
Il codice che si occupa della scrittura dell'header e dell'array di pixel è molto semplice:
Insieme di Mandelbrot
Cos'è un frattale appartenente all'insieme di Mandelbrot?
Esprimibile con la successione:
che ci dice se il numero complesso c appartiene al set (l'insieme di Mandelbrot) quando Zn tende all'infinito. Se appartiene al set viene tipicamente colorato di nero, altrimenti di bianco. Al bianco ci sono alternative, e lo si può colorare anche in base al numero di iterazioni (nell'immagine qui sopra viene privilegiato il blu e la colorazione dipende dalle iterazioni).
Nella formula Z e C appartengono all'insieme dei numeri complessi: essi sono composti da due parti, una reale ed una immaginaria. In rete si trovano numerose risorse a riguardo (e spiegate certamente meglio di quanto potrei fare io), mi limito solo a dire che consentono ad esempio di dare un valore alle radici negative; infatti
Ne segue che
Operazioni sui numeri complessi
Come detto nel paragrafo precedente, se abs(Zn) <= 2, allora C appartiene al set, e verrà quindi utilizzato il colore nero; diversamente verrà applicato un altro colore. Non sarà però sfuggito che le operazioni effettuate sono su numeri complessi e non su numeri reali.
La moltiplicazione tra due numeri complessi a e b avviene con questa formula:
La somma tra due numeri complessi a e b, avviene in questo modo:
Il valore assoluto di un numero invece lo si ottiene in questo modo (teorema di Pitagora):
a questo punto però possiamo evitare un calcolo dispendioso come la radice quadrata e confrontare abs(Z) non più con 2, bensì con 4.
Generare Frattali
Il frattale mostrato ad inizio articolo è stato generato dal codice qui sotto riportato.
Il cuore dell'algoritmo si compone di poco codice:
dopo aver ottenuto le coordinate "scalate", si procede con lo scoprire se C appartiene o no al set:
In base al valore di abs sappiamo quindi se C appartiene o no al set. A questo punto si sceglie il colore:
l'algoritmo in questo modo privilegia il colore blu sugli altri.
Tra le prime righe del file
in questo modo si può scegliere su quale parte andare a "zoomare". Il precedente frattale è stato generato con i valori riportati sopra (ma ritoccando un pò i colori, privilegiando verde e blu, ulteriormente).
Una zona interessante è la "seahorse valley", e la si può osservare settando questi valori:
Modificando i colori, privilegiando il verde invece del blu, e modificando anche il colore nero, si ottiene questo:
In una precedente versione un menu consentiva di impostare tutti i parametri; con l'ultimo aggiornamento ho apportato qualche modifica, e non ho ancora consentito una buona personalizzazione.
Ottimizzazioni: multithreading
Ciascun px di un frattale viene calcolato indipendentemente e questo rende piuttosto semplice suddividere l'immagine in parti più piccole e lasciare a più thread il compito di effettuare i calcoli.
Ho così implementato una funzione che consente di avere N quadranti, calcolando le coordinate di inizio e fine di ciascuno di essi.
Successivamente vengono avviati N thread; di default 4, ed il consiglio è di non superare i core fisici della macchina. Io ho un Notebook dotato di CPU Intel i7 6700HQ, con 4 core e 8 threads.
il punto di ingresso di ciascun thread è proprio la funzione calc() già analizzata in precedenza.
Con l'introduzione del multithreading i tempi di generazione si sono ridotti in modo considerevole:
Color Palette
La colorazione mostrata nel paragrafo precedente è in realtà incompleta; ho omesso una parte di codice:
E' possibile abilitare delle palette di colore; queste vengono caricate da un file di testo passato in input, composto semplicemente in questo modo:
la lunghezza è relativa alle righe del file, i colori sono espressi in decimale. Al momento della lettura del file questi valori vengono memorizzati in un 32bit (che è l'array chiamato palette); al momento della lettura del colore vengono quindi ottenute le 3 componenti richiamando get_channel().
Un esempio pratico può essere questo:
Si tratta dei colori che compongono l'immagine mostrata da Wikipedia. Decommentando la direttiva PALETTE in mandelbrot.h, si abilita questa funzionalità.
Il risultato è mostrato qui sotto:
Conclusione
Spero sia stato un viaggio breve ma interessante!
Il "progetto" lo trovate sul mio GitHub, sotto il nome di fraCtal.
Ovviamente qualora voleste scrivere una vostra versione in qualsiasi linguaggio e pubblicarla in questo topic, assieme magari a qualche immagine, sarebbe ben accetta. :)
Questa volta voglio parlare di qualcosa di insolito rispetto al mio target. Anni fa scrissi un piccolissimo generatore di frattali, che pubblicai più avanti su Github, appartenendi all'insieme di Mandelbrot. Nulla di particolarmente avanzato e senza l'utilizzo di tecniche particolari di colorazione, sia chiaro; circa 1 mese fa ho aggiunto il supporto al multithreading.
Ho deciso di scrivere un brevissimo articolo solo per incuriosire chi magari non programma da molto tempo e non conosce i frattali. Spero anche di suscitare l'interese dei i più navigati, ovviamente!
Si trattava di una partecipazione ad un topic su un forum quella che mi fece scrivere queste due righe, quindi decisi come altri di avere poche dipendenze. Le immagini vengono generate in PPM; quelle che vedete nell'articolo sono quindi state tutte convertite.
L'articolo avrà il seguente svolgimento:
- Immagini PPM: cosa sono?
- Insieme di Mandelbrot
- Operazioni sui numeri complessi
- Generare Frattali
- Ottimizzazioni: multithreading
- Color Palette
Immagini PPM: cosa sono?
Portable pixmap format, conosciuto anche come PPM, è un formato di immagine definito nel progetto Netpbm: allo stesso progetto appartengono PGM e PBM rispettivamente immagini in scala di grigi e in bianco e nero. PPM invece permette di rappresentare immagini a colori (0-255) in formato RGB. Sono molto semplici da gestire e possono essere in ASCII o binarie.
Codice:
P6
# commento
w h
255quello riportato qui sopra è l'header di un'immagine PPM rappresentata in binario (P6) larga w pixel e alta h pixel; ogni pixel dell'immagine ha una profondità di 24bit (8 * 3, dove 3 sono i canali di colore, RGB).
Il codice che si occupa della scrittura dell'header e dell'array di pixel è molto semplice:
C:
void save_image(char *imageName, char* pixels, size_t w, size_t h)
{
FILE *file = fopen(imageName, "wb");
fprintf(file, "%s%d%s%d%s", "P6\n#DispatchCode PPMImageLib\n", w, " ", h, "\n255\n");
fwrite(pixels, sizeof(char), 3+(w*h*3), file);
fclose(file);
}Insieme di Mandelbrot
Cos'è un frattale appartenente all'insieme di Mandelbrot?
Esprimibile con la successione:
che ci dice se il numero complesso c appartiene al set (l'insieme di Mandelbrot) quando Zn tende all'infinito. Se appartiene al set viene tipicamente colorato di nero, altrimenti di bianco. Al bianco ci sono alternative, e lo si può colorare anche in base al numero di iterazioni (nell'immagine qui sopra viene privilegiato il blu e la colorazione dipende dalle iterazioni).
Nella formula Z e C appartengono all'insieme dei numeri complessi: essi sono composti da due parti, una reale ed una immaginaria. In rete si trovano numerose risorse a riguardo (e spiegate certamente meglio di quanto potrei fare io), mi limito solo a dire che consentono ad esempio di dare un valore alle radici negative; infatti
sqrt(-1) = i, vale a dire l'unità immaginaria, esprimibile anche come la coppia ordinata (0, 1i). In questo caso 0 è la componente reale, mentre 1i è la componente immaginaria del numero.Ne segue che
Z = a + ib. A questo punto occorre però limitare la serie se vogliamo scrivere un algoritmo; ecco che sempre la matematica ci viene in soccorso: fissato un numero N di iterazioni, possiamo dire che se il valore assoluto di Z è > di 2, allora non tornerà più inferiore a 2, ed il loop si può interrompere; in tal caso verrà applicato un colore in quanto c non apparterrà all'insieme.Operazioni sui numeri complessi
Come detto nel paragrafo precedente, se abs(Zn) <= 2, allora C appartiene al set, e verrà quindi utilizzato il colore nero; diversamente verrà applicato un altro colore. Non sarà però sfuggito che le operazioni effettuate sono su numeri complessi e non su numeri reali.
La moltiplicazione tra due numeri complessi a e b avviene con questa formula:
C:
re = a.re * b.re - a.im * b.im;
im = a.re * b.im + a.im * b.re;La somma tra due numeri complessi a e b, avviene in questo modo:
C:
re = a.re + b.re;
im = a.im + b.im;Il valore assoluto di un numero invece lo si ottiene in questo modo (teorema di Pitagora):
C:
abs = sqrt(a.re * a.re + a.im * a.im)a questo punto però possiamo evitare un calcolo dispendioso come la radice quadrata e confrontare abs(Z) non più con 2, bensì con 4.
Generare Frattali
Il frattale mostrato ad inizio articolo è stato generato dal codice qui sotto riportato.
Il cuore dell'algoritmo si compone di poco codice:
C:
void *calc(_square_coords *coords)
{
for(int i = coords->start_y; i < coords->end_y; i++)
{
double scaled_y = Y1 + (((double)i * DIFF_Y) / (double)coords->h);
for(int j = coords->start_x; j<coords->end_x; j++)
{
double scaled_x = X1 + (((double)j * DIFF_X) / (double)coords->w);
complex z, c;
z.re = scaled_x;
z.im = scaled_y;
c.re = scaled_x;
c.im = scaled_y;
int k;
double abs=0.0;
for(k = 0; k < coords->max_iterations; k++)
{
z = complex_sum(complex_mul(z, z), c);
abs = complex_abs(z);
if(abs > 4.0) {
break;
}
}
int index = linear_index(j, i, coords->w);
if(abs <= 4.0 )
{
pixels[index] = 0;
pixels[index+1] = 0;
pixels[index+2] = 0;
}
else
{
int r = k >> 2;
int g = k >> 1;
int b = k << 3;
r = (r > 128) ? 128 : r;
g = (g > 128) ? 128 : g;
b = (b > 255) ? 255 : b;
pixels[index ] = r;
pixels[index + 1] = g;
pixels[index + 2] = b;
}
}
}
}
C:
double scaled_y = Y1 + (((double)i * DIFF_Y) / (double)coords->h);
double scaled_x = X1 + (((double)j * DIFF_X) / (double)coords->w);dopo aver ottenuto le coordinate "scalate", si procede con lo scoprire se C appartiene o no al set:
C:
complex z, c;
z.re = scaled_x;
z.im = scaled_y;
c.re = scaled_x;
c.im = scaled_y;
int k;
double abs=0.0;
for(k = 0; k < coords->max_iterations; k++)
{
z = complex_sum(complex_mul(z, z), c);
abs = complex_abs(z);
if(abs > 4.0) {
break;
}
}In base al valore di abs sappiamo quindi se C appartiene o no al set. A questo punto si sceglie il colore:
C:
int index = linear_index(j, i, coords->w);
if(abs <= 4.0 )
{
pixels[index] = 0;
pixels[index+1] = 0;
pixels[index+2] = 0;
}
else
{
int r = k >> 2;
int g = k >> 1;
int b = k << 3;
r = (r > 128) ? 128 : r;
g = (g > 128) ? 128 : g;
b = (b > 255) ? 255 : b;
pixels[index ] = r;
pixels[index + 1] = g;
pixels[index + 2] = b;
}l'algoritmo in questo modo privilegia il colore blu sugli altri.
Tra le prime righe del file
mandelbrot.c del mio piccolo progetto compaiono alcune direttive:
C:
#define X0 (2.0)
#define Y0 (1.7)
#define X1 (-2.0)
#define Y1 (-1.7)
#define DIFF_Y (Y0-Y1)
#define DIFF_X (X0-X1)in questo modo si può scegliere su quale parte andare a "zoomare". Il precedente frattale è stato generato con i valori riportati sopra (ma ritoccando un pò i colori, privilegiando verde e blu, ulteriormente).
Una zona interessante è la "seahorse valley", e la si può osservare settando questi valori:
C:
#define X0 (0.37)
#define Y0 (0.21)
#define X1 (0.40)
#define Y1 (0.26)
Modificando i colori, privilegiando il verde invece del blu, e modificando anche il colore nero, si ottiene questo:
In una precedente versione un menu consentiva di impostare tutti i parametri; con l'ultimo aggiornamento ho apportato qualche modifica, e non ho ancora consentito una buona personalizzazione.
Ottimizzazioni: multithreading
Ciascun px di un frattale viene calcolato indipendentemente e questo rende piuttosto semplice suddividere l'immagine in parti più piccole e lasciare a più thread il compito di effettuare i calcoli.
Ho così implementato una funzione che consente di avere N quadranti, calcolando le coordinate di inizio e fine di ciascuno di essi.
Successivamente vengono avviati N thread; di default 4, ed il consiglio è di non superare i core fisici della macchina. Io ho un Notebook dotato di CPU Intel i7 6700HQ, con 4 core e 8 threads.
C:
void generation(mandelbrot_info mandelbrot, char *filename)
{
pixels = (char*) calloc(3 + (mandelbrot.w * mandelbrot.h * 3), 8);
pthread_t square[mandelbrot.n_quadrants];
for(int i=0; i<mandelbrot.n_quadrants; i++)
pthread_create(&square[i], NULL, (void*) calc, (void*)&mandelbrot.coords[i]);
for(int i=0; i<mandelbrot.n_quadrants; i++)
pthread_join(square[i], NULL);
free(mandelbrot.coords);
save_image(filename, pixels, mandelbrot.w, mandelbrot.h);
free(pixels);
}il punto di ingresso di ciascun thread è proprio la funzione calc() già analizzata in precedenza.
Con l'introduzione del multithreading i tempi di generazione si sono ridotti in modo considerevole:
| 1-thread | 4-thread | (x0, x1, y0, y1) |
| 166.30s | 57.37s | (0.37, 0.40, 0.21, 0.16) |
| 198.322s | 111.56s | (0.37, 0.40, 0.21, 0.26) |
| 46.11s | 21.60s | (2.0, -2.0, 1.7, -1.7) |
Color Palette
La colorazione mostrata nel paragrafo precedente è in realtà incompleta; ho omesso una parte di codice:
C:
#ifndef PALETTE
int r = k >> 2;
int g = k >> 1;
int b = k << 3;
r = (r > 128) ? 128 : r;
g = (g > 128) ? 128 : g;
b = (b > 255) ? 255 : b;
pixels[index ] = r;
pixels[index + 1] = g;
pixels[index + 2] = b;
#else
uint32_t color = palette[k % palette_size];
pixels[index ] = get_channel(color,16);
pixels[index + 1] = get_channel(color,8);
pixels[index + 2] = get_channel(color,0);
#endifE' possibile abilitare delle palette di colore; queste vengono caricate da un file di testo passato in input, composto semplicemente in questo modo:
Codice:
length
r,g,b
r,g,b
...la lunghezza è relativa alle righe del file, i colori sono espressi in decimale. Al momento della lettura del file questi valori vengono memorizzati in un 32bit (che è l'array chiamato palette); al momento della lettura del colore vengono quindi ottenute le 3 componenti richiamando get_channel().
Un esempio pratico può essere questo:
Codice:
16
66,30,15
25,7,26
9,1,47
4,4,73
0,7,100
12,44,138
24,82,177
57, 125,209
134,181,229
211,236,248
241,233,191
248,201,95
255,170,0
204,128,0
153,87,0
106,52,3Si tratta dei colori che compongono l'immagine mostrata da Wikipedia. Decommentando la direttiva PALETTE in mandelbrot.h, si abilita questa funzionalità.
Il risultato è mostrato qui sotto:
Conclusione
Spero sia stato un viaggio breve ma interessante!
Il "progetto" lo trovate sul mio GitHub, sotto il nome di fraCtal.
Ovviamente qualora voleste scrivere una vostra versione in qualsiasi linguaggio e pubblicarla in questo topic, assieme magari a qualche immagine, sarebbe ben accetta. :)
