Discussione Disegnare frattali: l'insieme di Mandelbrot

[DispatchCode]

No, niente esercizi. È più un divertimento, questa volta (per chi vorrà partecipare)

Ho visto che se ne è parlato in passato, anni fa, ma volevo rispolverare l'argomento visti anche i nuovi utenti, specie i più giovani e "nuovi" della programmazione, che magari non conoscono ancora l'argomento.

Si tratta dei frattali di Mandelbrot, "creature" come questa:

Potete anche eseguire lo "zoom" nel frattale, e vedere sezioni come questa:

La si può anche rendere più definita:

Trovate tutte le info qui: https://en.wikipedia.org/wiki/Mandelbrot_set

Non fatevi spaventare dalle formule, si tratta alla fine di un paio di cicli for innestati e alcune operazioni sui numeri complessi. Su Wiki potete trovare lo pseudocodice per aiutarvi con l'implementazione.

Potete implementare l'algoritmo con qualsiasi linguaggio vogliate; io utilizzai C. Non ho una GUI, magari in futuro la implementerò. Genero un'immagine PNG utilizzando lodepng.

Attendo a condividere il sorgente per vedere prima qualche vostro intervento.

 

[DispatchCode]

Visto che qualche giorno è trascorso, pubblico il codice. Ho rimosso alcune parti futili:

void *calc(_square_coords *coords)
{
    double x_step = DIFF_X / coords->w;
    double y_step = DIFF_Y / coords->h;

    // Blue and red from 0.00 to 255.0
    double blue_a = 0.0, blue_b = 255.0;
    double red_a  = 0.0, red_b  = 255.0;

    // Fractal must be red on x axis and blue on y axis
    double rstepx = (red_b - red_a) / coords->w;
    double bstepy = (blue_b - blue_a) / coords->h;

    for(int i = coords->start_y; i < coords->end_y; i++)
    {
        complex z, c;
        c.im = Y1 + (((double)i) * y_step);

        for(int j = coords->start_x; j<coords->end_x; j++)
        {
            c.re = X1 + (((double)j) * x_step);
            z.re = 0;
            z.im = 0;

            int k = 0;
            double abs = 0.0;
            while(k++ < coords->max_iterations && abs < 4.0)
            {
                z   = complex_sum(complex_mul(z, z), c);
                abs = complex_abs(z);
            }

            int index = linear_index(j, i, coords->w);

            int r = (int)(red_a + rstepx * j);
            int b = (int)(blue_a + bstepy * i);

            int val_in = 0, val_out = 0;
            int g = 0;
            if(k >= 32) {
                g = (int)(255.0 - ((((double)k + val_in)) / (coords->max_iterations + val_in) * (255.0-val_out)));
            }

            double atten = ((double)k / coords->max_iterations);
            atten = 1.0 - (atten*atten);

            r = (int)((double)r * atten);
            b = (int)((double)b * atten);

            pixels[index    ] = (char)r;
            pixels[index + 1] = (char)g;
            pixels[index + 2] = (char)b;
            pixels[index + 3] = (char)255; // alpha channel

        }
    }
}

Questo è il cuore dell'algoritmo.
Il codice si compone di altre parti in realtà:

void init(size_t w, size_t h, size_t n_threads,size_t iterations, mandelbrot_info *info)
{
    size_t middle_x = w >> 1;
    size_t middle_y = h >> 1;

    info->h = h;
    info->w = w;
    info->n_quadrants = n_threads;
    info->coords      = calloc(n_threads, sizeof(_square_coords));

    int m_th = n_threads >> 1;
    for(int i=0; i < n_threads; i++)
    {
        if(i < m_th)
        {
            info->coords[i].start_x = 0;
            info->coords[i].end_x = middle_x;
        }
        else
        {
            info->coords[i].start_x = middle_x;
            info->coords[i].end_x = w;
        }

        if(i&1)
        {
            info->coords[i].start_y = middle_y;
            info->coords[i].end_y = h;
        }
        else
        {
            info->coords[i].start_y = 0;
            info->coords[i].end_y = middle_y;
        }

        info->coords[i].w = w;
        info->coords[i].h = h;
        info->coords[i].max_iterations = iterations;
    }

#ifdef PALETTE
   load_palette("palette.txt");
#endif
}

/*
 * Run threads and wait until they ended
 *
 */
void generation(mandelbrot_info mandelbrot, char *filename)
{
    pixels = (char*) calloc(4 + (mandelbrot.w * mandelbrot.h * 4), 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 png image using lodepng
    unsigned error = lodepng_encode32_file(filename, pixels, mandelbrot.w, mandelbrot.h);
    if(error) {
        printf("error %u: %s\n", error, lodepng_error_text(error));
    }

    free(pixels);
}

Avevo ottimizzato il codice utilizzando il multi-threading: l'immagine viene suddivisa in 4 quadranti, e ciascuno viene calcolato separatamente (il valore, 4, può essere cambiato dal main).

Il mio codice completo con tanto di spiegazione e svariate immagini lo trovate qui: https://github.com/DispatchCode/fraCtal

 

[CrazyMonk]

Mi scuso per il ritardo, purtroppo sono sempre impegnato e non ho molto tempo libero.
Io ho provato a fare qualcosa del genere utilizzando Python e le librerie numpy e matplotlib:

import numpy
import matplotlib.pyplot as plt

#Creo una matrice nella quale si stabilisce una range di coordinate
#Utilizzando la funzione linspace() genero un numero di punti equidistanti per righe e colonne
#Restituisce dei vettori bidimensionali
def crea_matrice(x_min, x_max, y_min, y_max, densita_pixel):
    riga = numpy.linspace(x_min, x_max, int((x_max - x_min) * densita_pixel))
    colonna = numpy.linspace(y_min, y_max, int((y_max - y_min) * densita_pixel))
    return riga[numpy.newaxis, :] + colonna[:, numpy.newaxis] *1j

#Verifico se ogni numero generato appartiene al set di Mandelbrot attraverso la semplice regola del raggio minore/uguale di 2
#Restituisco una sequenza di true/false
def appartieneSetMandelbrot(c, iterazioni):
    n = 0
    for i in range(iterazioni):
        n = ((n**2) + c)
    return abs(n) <= 2

#Attraverso la sequenza di true/false trovata al passo precedente posso restituire un array monodimensionale contenente solo i numeri che appartengono al set di mandelbrot
def filtraMembriAppartenenti(c, iterazioni):
    mask = appartieneSetMandelbrot(c, iterazioni)
    return c[mask]

#Definisco i parametri per il disegno del frattale
densita_pixel = 512
iterazioni = 20
matrice = crea_matrice(-2, 0.5, -1.5, 1.5, densita_pixel)
plt.imshow(appartieneSetMandelbrot(matrice, iterazioni), cmap="binary")
plt.gca().set_aspect("equal")
plt.axis("off")
plt.tight_layout()
plt.show()

P.s : Una cosa che ho capito, però, è che la generazione dei numeri di mandelbrot cresce esponenzialmente anche con "c" molto piccoli, ma non ho capito come limitare questa crescita esponenziale per evitare errori di overflow.

EDIT: Questo è il risultato grafico

 

[DispatchCode]

Figurati, anzi, grazie per aver partecipato!

L'ho provato anche io, mi suona strano quell'overflow in effetti. Puoi provare ad effettuare il controllo su abs(n) ad ogni iterazione? Come condizione del ciclo insomma.
Come ho fatto io qui praticamente: while(k++ < coords->max_iterations && abs < 4.0).

Prova a verificare il tipo di dato che stai usando con numpy: probabilmente puoi impostare un tipo di dato più grande quando crei la matrice (32/64bit).

Bel lavoro comunque, rimango sempre stupito dalla velocità con cui si fanno le cose in Python. xD

Ma la funzione filtraMembriAppartenenti() è un refuso, giusto? Mi sembra la si possa togliere senza problemi.

PS.
Ah aggiungo un piccolo HINT: se invece di restituire true/false restituisci il numero di iterazioni, puoi anche "sfumare" in base al valore che ottieni. Una volta usavo questo algoritmo: https://github.com/DispatchCode/fra...d66158c41bb0f6e293115e142e5c18c555d656e47R128

Anche nella versione attuale uso k per colorare, in realtà. Per l'algoritmo attuale mi aveva dato prezioni hint BrutPitt (è un gran appassionato, se gli guardi il profilo github vedrai un progetto... che fa paura; ed è un fisico che fa l'informatico da anni, quindi con le formule si trova a proprio agio).

 

[CrazyMonk]

Figurati, anzi, grazie per aver partecipato!

L'ho provato anche io, mi suona strano quell'overflow in effetti. Puoi provare ad effettuare il controllo su abs(n) ad ogni iterazione? Come condizione del ciclo insomma.
Come ho fatto io qui praticamente: while(k++ < coords->max_iterations && abs < 4.0).

Prova a verificare il tipo di dato che stai usando con numpy: probabilmente puoi impostare un tipo di dato più grande quando crei la matrice (32/64bit).

Bel lavoro comunque, rimango sempre stupito dalla velocità con cui si fanno le cose in Python. xD

Ma la funzione filtraMembriAppartenenti() è un refuso, giusto? Mi sembra la si possa togliere senza problemi.

PS.
Ah aggiungo un piccolo HINT: se invece di restituire true/false restituisci il numero di iterazioni, puoi anche "sfumare" in base al valore che ottieni. Una volta usavo questo algoritmo: https://github.com/DispatchCode/fra...d66158c41bb0f6e293115e142e5c18c555d656e47R128

Anche nella versione attuale uso k per colorare, in realtà. Per l'algoritmo attuale mi aveva dato prezioni hint BrutPitt (è un gran appassionato, se gli guardi il profilo github vedrai un progetto... che fa paura; ed è un fisico che fa l'informatico da anni, quindi con le formule si trova a proprio agio).

Si dovrei provare a fare un controllo sul ciclo, hai ragione. In effetti l'overflow è abbastanza giustificato perché i numeri, anche con valori di c piccoli, crescono a dismisura. Questi, per esempio, sono i valori che saltano fuori con c = 1, iterando 10 volte:

1
2
5
26
677
458330
210066388901
44127887745906175987802
1947270476915296449559703445493848930452791205
3791862310265926082868235028027893277370233152247388584761734150717768254410341175325352026

Poi ci sono valori particolari per i quali il valore dell'n-esimo numero oscilla sempre tra due estremi, per esempio c = -1, iterando sempre 10 volte:

-1
0
-1
0
-1
0
-1
0
-1
0

Per quanto riguarda la funzione filtraMembriAppartenenti() è un refuso, infatti, non serve. Il python ti permette di fare cose complesse con, relativamente, poche righe di codice. E' eccezionale...soprattutto per la vastità di librerie che possiede, in qualsiasi campo. Molte librerie utilizzano la vettorizzazione, per esempio, quindi puoi fare cose assai grandi con poco codice: il problema poi diventa la capacità di astrazione però ahahaha. E' un po' come il matlab, non so se conosci questo linguaggio, con cui puoi fare operazioni su matrici multidimensionali con veramente pochi comandi ed eliminando l'utilizzo dei cicli: è un casino, però, sviluppare la logica applicativa perché ti allontani molto dall'approccio della programmazione procedurale, la quale ragiona su singole componenti, poi il codice diventa anche poco leggibile ecc.

Grazie anche del consiglio sulle sfumature, proverò ad applicarlo. Sono molto curioso di andare a visitare la pagina di questo BrutPitt...i fisici sono pieni di sorprese

 

[DispatchCode]

Si, ho presente Matlab. Mai utilizzato però, il mio campo di applicazione è un altro, quindi non mi è mai servito.
Con Python avevo smanettato un pò anni fa, c'era ancora la versione 2, ma niente di serio.

Prova ad aggiungere quella condizione, magari riesci già a rimediare al problema overflow. Inoltre puoi anche limitare il ciclo ad abs < 4 invece di 2.

In merito a BrutPitt (Michele) il suo progetto è questo: https://github.com/BrutPitt/glChAoS.P
Buon divertimento xD

 

[pvssygino]

Ciao! Ho notato che i frattali sono spesso considerati un argomento strettamente matematico e astratto, a causa della presenza di numeri complessi. Tuttavia, come ingegnere elettronico, ti posso assicurare che questi oggetti sono utilizzati anche in applicazioni pratiche, come la creazione di antenne o il miglioramento della densità dei collegamenti su schede PCB. In particolare, le antenne frattali sono costituite da elementi ripetitivi di dimensioni diverse, che consentono di migliorare le prestazioni dell'antenna stessa. Al contrario di un singolo filo, gli elementi frattali sono in grado di catturare e trasmettere segnali elettromagnetici con maggiore efficienza e precisione. Ecco una foto che mostra un'antenna frattale

Son sicuro che avrete visto in giro un'antenna del genere, spesso i traingolo alla base hanno un perimetro maggiore di quelli alla punta. Il triangolo utilizzato è quello di Sierpiński, in questo caso con iterazione = 2, potete trovare questo frattale nella documentazione di Wolfram ( Spierpniski ).

 

[CrazyMonk]

Ciao! Ho notato che i frattali sono spesso considerati un argomento strettamente matematico e astratto, a causa della presenza di numeri complessi. Tuttavia, come ingegnere elettronico, ti posso assicurare che questi oggetti sono utilizzati anche in applicazioni pratiche, come la creazione di antenne o il miglioramento della densità dei collegamenti su schede PCB. In particolare, le antenne frattali sono costituite da elementi ripetitivi di dimensioni diverse, che consentono di migliorare le prestazioni dell'antenna stessa. Al contrario di un singolo filo, gli elementi frattali sono in grado di catturare e trasmettere segnali elettromagnetici con maggiore efficienza e precisione. Ecco una foto che mostra un'antenna frattale

Son sicuro che avrete visto in giro un'antenna del genere, spesso i traingolo alla base hanno un perimetro maggiore di quelli alla punta. Il triangolo utilizzato è quello di Sierpiński, in questo caso con iterazione = 2, potete trovare questo frattale nella documentazione di Wolfram ( Spierpniski ).

Forte! Non credevo che i frattali avessero anche delle applicazioni pratiche, quindi sono anche utili, oltre ad essere molto interessanti dal punto di vista artistico/matematico.