The Xna-Way: Tutorial 9: Fractal Generator

Come avevo preannunciato nell'ultimo post, ho realizzato un piccolo componente che genera alcuni frattali.
Il tutto è scritto in .Net ed utilizza gli strumenti forniti da XNA per creare e manipolare le immagini usate per la generazione dei frattali.
Il "programma" che illustrerò è abbastanza lunghetto rispetto agli altri, ma vedrò nel possibile di spiegarlo al meglio delle mie capacità.



Come sempre vi esorto a farmi notare se ho detto delle boiate immani!

Nota: nell'ultimo post avevo fatto vedere diverse immagini di frattali realizzate con il codice che andrò ad illustrare, però il codice che vi fornirò adesso potrà generare un solo tipo di frattale, e cioè il frattale di Mandelbrot.
Gli altri non li metterò per ora, anche se il codice è molto simile a quello che userò qua.

Intanto cominciamo col dire come ho pensato il tutto:
>il calcolo che dobbiamo fare è dato dalla successione:
z(n+1) = z(n)^2 + c
dove z è il valore del punto che vogliamo calcolare, e c è una costante che dobbiamo sommare di volta in volta
per z(0) valore sarà c
per z(1) il valore sarà c^2 + c
e così via.

>i punti su cui calcoliamo la successione sono quelli del piano complesso, quindi dovremo fare tutti i calcoli utilizzando l'algebra dei numeri complessi

>il valore di c cambia a seconda di quale punto del piano complesso stiamo calcolando il valore

>un punto appertiene all'insieme di Mandelbrot se Zn è in modulo minore di 2, dove n indica per quante volte dobbiamo ripetere il calcolo prima di considerare il punto in questione appartenente all'insieme. Quindi se dopo un numero i di ripetizioni del calcolo della successione il valore supera 2, il punto non appartiene all'insieme e verrà colorato in modo particolare, mentre se il punto appartiene all'insieme sarà colorato di nero.

NOTA: per la gestione dei calcoli con i numeri complessi ho usato una classe che ho trovato su internet, del quale ho fatto brutalmente copy-paste. Lo che di solito non si fa così, ma non era la parte dui numeri complessi che mi interessava implementare. Sicuramente l'uso di questa classe può aggiungere un po' di tempo di calcolo in più, quindi per velocizzare avrei potuto implementare il calcolo a manina, ma la "complessità" dei numeri complessi (scusate il gioco di parole XD) mi ha tolto la voglia.

Avrò quindi bisono di sapere in quale zona del piano complesso sto andando a calcolare la successione, quindi dovrò definire il limite inferiore e superiore per l'asse reale (considerato l'asse delle ascisse) e l'asse complesso (considerato delle ordinate).
Poi dovrò sapere quanto è lungo il lato delle ascisse e delle ordinate del piano (cioè la differenza tra il limite superiore ed inferiore degli assi), dato che questi valori poi mi servono per calcolare il valore di incremento. Si deve considerare che la texture (che rappresenta il frattale) sarà una "approssimazione" del piano complesso, e dato che non posso avere tanti pixel quanti sono i punti del piano (dato che sarebbero infiniti) devo assegnare ad ogni pixel un punto del piano e calcolare il valore per quel punto, e per fare ciò, quando passo da un pixel al successivo, devo sapere di quando avando/devo incrementare per passare al punto successivo.

Spieghiamo meglio questa cosa dei pixel...
Come forse saprete i punti del piano sono infiniti, mentre per quanto possano essere piccoli e fitti i pixel di un'immagine essi sono in numero ben finito.
Se ipotizzo i pixel dell'immagine come una "discretizzazione" del piano complesso (passatemi questo brutto termine che a molti non piace), allora devo assegnare ad ogni pixel un certo punto del piano, fare i calcoli su quel punto, e colorare il pixel a seconda del risultato ottenuto.
Ora prendiamo il primo pixel dell'immagine, quello che sta in alto a sinistra per interci (perchè nelle immagini l'asse delle X va da sinistra a destra e quello delle Y dall'alto verso il basso), e diamogli come punto quello identificato dall'estremo inferiore per l'asse X e dall'estremo superiore per l'alle Y. Quando poi vogliamo passare al punto successivo (per ora considero lo spostamento solo sull'asse reale, cioè delle X) dovrò incrementare di una certa quantità.
Questa quantità che devo sommare per passare da un punto all'altro la posso calcolare sapendo quanto è lungo il lato reale che sto prendendo in cosiderazione (quindi estremo superiore di X - estremo inferirio di X) e va diviso per la larghezza dell'immagine (in pixel naturalmente). In questo modo per passare da un punto a l'altro si fa una semplice somma.
Notate che la stessa cosa va fatta/calcolta per l'asse delle Y, trovando la quantità di incremento dividendo il lato verticale del piano per il numero di pixel dell'altezza.

Poi ho bisongno del massimo numero di iterazioni da calcolare prima di considerare il punto come appartenente all'insieme. Ricordo che un punto viene considerato parte dell'insieme di Madelbrot e quindi colorato di nero, solo se dopo N iterazioni (con N numero massimo di iterazioni) del calcolo della successione, il valore rispetta ancora la legge.
Naturalmente avrò bisogno anche di altre variabili: per esempio dovrò poter memorizzare il valore di ogni punto del piano di cui sto calcolando la successione, il numero di iterazioni fatte per ogni punto calcolato, e intanto che ci siamo anche per quali indici/pixel/punto devo ancora continuare a calcolare la successione.
Le variabili definite a questo scopo sono:

  protected double infx = -2, supx = 2, infy = -2, supy = 2; //limiti del piano complesso su cui calcoare la successione
  protected double latox, latoy, incx, incy; //lunghezza dei lati e di quando devo incrementare il valore 
          //da un pixel e l'altro
  protected int maxIteration = 250; //numero massimo di iterazione

  protected Texture2D fractal; //texture
  protected bool complete = true; //indica se il calcolo è completo o no
  protected int width, height; //altezza e larghezza della texture
  protected Color[] pixData; //vettore contenente il colore dei pixel della texture
  protected int total; // numero totale dei pixel della texture

  protected Queue<int> indici = new Queue<int>(); //indice del pixel che devo calcolare
  protected int[,] iterazioni; //numero di iterazioni per ogni punto
  protected Cmplx[,] valori; //valore di ogni punto
  protected Cmplx zCm = new Cmplx(), cCm = new Cmplx(); //numeri complessi per il calcolo

NOTA: il perchè ho usato una queue (struttura dati che implementa una coda, cioè una "lista" di elementi con estrazione in testa testa ed inserimento in coda) sarà chiaro in seguito.

La classe MandelBrotFractal è figlia di GameComponent, in questo modo potrà avere il suo metodo update e eseguire parte del calcolo ad ogni chiamata del metodo. Inoltre implementerà un'interfaccia da me definita che esporrà i metodi per agire sul frattale.

public interface FractalI
{
  Texture2D Fractal { get; }
  void resizeTexture(int width, int height);
  bool Complete { get; }
  int Width { get; }
  int Heigth { get; }
  void UpdateFractalDouble(double infx, double supx, double infy, double supy, int MAX);
}

Il significato mi pare abbastanza esplicito, comunque quando si vedrà come vengono usati i metodi sarà ancora più chiaro.

Il codice completo della classe che genera il MandelBrot è questo:
Show/Hide

public class MandelBrotFractal : Microsoft.Xna.Framework.GameComponent, FractalI
{
    protected double infx = -2, supx = 2, infy = -2, supy = 2; //limiti del piano complesso su cui calcoare la successione
    protected double latox, latoy, incx, incy; //lunghezza dei lati e di quando devo incrementare il valore 
            //da un pixel e l'altro
    protected int maxIteration = 250; //numero massimo di iterazione

    protected Texture2D fractal; //texture
    protected bool complete = true; //indica se il calcolo è completo o no
    protected int width, height; //altezza e larghezza della texture
    protected Color[] pixData; //vettore contenente il colore dei pixel della texture
    protected int total; // numero totale dei pixel della texture

    protected Queue<int> indici = new Queue<int>(); //indice del pixel che devo calcolare
    protected int[,] iterazioni; //numero di iterazioni per ogni punto
    protected Cmplx[,] valori; //valore di ogni punto
    protected Cmplx zCm = new Cmplx(), cCm = new Cmplx(); //numeri complessi per il calcolo

    public MandelBrotFractal(Game game, double infx, double supx, double infy, double supy, int maxInt)
        : base(game)
    {
        //prendo i valori della larghezza e altezza della texture
        width = this.Game.GraphicsDevice.PresentationParameters.BackBufferWidth;
        height = this.Game.GraphicsDevice.PresentationParameters.BackBufferHeight;
        //numero totale di pixel
        total = width * height;
        pixData = new Color[total]; //creo l'array per i colori dei pixel

        //creo la texture
        fractal = new Texture2D(this.Game.GraphicsDevice,
            width,
            height, 1,
            TextureUsage.None, SurfaceFormat.Color);
        //creo la matrice che conterrà il numero di iterazioni fatte per ogni punto
        iterazioni = new int[height, width];
        //creo la matrice di numeri complessi che conterrà il valore di ogni punto del piano che sto calcolando
        //e di volta in volta il calcolo parziale
        valori = new Cmplx[height, width];
        for (int i = 0; i < height; i++)
            for (int j = 0; j < width; j++)
                valori[i, j] = new Cmplx();
    }

    public override void Initialize()
    {
        base.Initialize();
    }

    public override void Update(GameTime gameTime)
    {
        if (!complete)
        {
            //se il frattale è da completare eseguo il calcolo
            calcolaFrattale();
        }
        base.Update(gameTime);
    }

    protected virtual Cmplx calcola(Cmplx z, Cmplx c)
    {
        //calcolo il valore
        Cmplx ris = z.sqr().add(c);
        return ris;
    }

    protected Cmplx ris;

    private void calcolaFrattale()
    {
        //variabili usate per inizializzare c
        double cre, cim;

        int indice = 0;
        //memorizzo quanti indici devo calcolare
        int count = indici.Count;
        for(int i = 0; i < count; i++) //faccio count estrazioni dalla coda
        {
            indice = indici.Dequeue();
            //calcolo il corrispettivo indice di riga e di colonna
            int r = indice / width;
            int c = indice % width;
            //calcolo il valore di c da usare (potrei memorizzarlo in una matrice?)
            cre = infx + (c * incx);
            cim = infy + (r * incy);

            //recupero il valore di z
            zCm = valori[r, c];
            //controllo che il valore sia in modulo minore di 4 e che non abbia già fatto il max num di iterazioni
            if (zCm.modsq() <= 4 && iterazioni[r, c] < maxIteration)
            {
                cCm.set(cre, cim);
                ris = calcola(zCm, cCm); //eseguo il calcolo
                valori[r, c] = ris; //aggiorno
                iterazioni[r, c]++; //aggiorno
                indici.Enqueue(indice); //reinserisco l'indice nella coda per il prossimo calcolo
            }
            else
            {
                if (iterazioni[r, c] == maxIteration)
                {
                    //il punto fa parte dell'insieme: lo coloro di nero
                    pixData[indice] = Color.Black;
                }
                else
                {
                    //il punto non fa parte dell'insieme
                    float val = (float)iterazioni[r, c] / maxIteration;
                    int h, s, v;
                    HsvToRgb(360.0f * val, 1.0f, 1.0f, out h, out s, out v);
                    //lo coloro in modo opportuno
                    pixData[indice] = new Color((Byte)h, (Byte)s, (Byte)v);

                    //Nota: cambiando come si calcola il colore posso ottenere scale di colore diverse
                }
            }
        }
        try
        {
            fractal.SetData(pixData);//aggiorno la texture
        }
        catch (Exception e)
        { }
        if (indici.Count == 0)
        {
            //se la coda risulta vuota vuol dire che ho finito di calcolare
            this.Game.Window.Title = "Render complete";
            complete = true;
        }
    }

    private Color makeColor(float r, float g, float b)
    {
        System.Drawing.Color color = System.Drawing.Color.FromArgb((int)(r), (int)(g), (int)(b));
        return new Color(color.R, color.G, color.B);
    }

//qui c'è la funzione per il passaggio dal HSV a RGB
//guardate i sorgenti

    int Clamp(int i)
    {
        if (i < 0) return 0;
        if (i > 255) return 255;
        return i;
    }

    #region FracalI Membri di

    public Texture2D Fractal
    {
        get { return fractal; }
    }

    /// <summary>
    /// aggiorna il calcolo del frattale dando gli estremi del piano su cui calcolare e il massimo numero di iterazioni
    /// </summary>
    /// <param name="infx"></param>
    /// <param name="supx"></param>
    /// <param name="infy"></param>
    /// <param name="supy"></param>
    /// <param name="maxIt"></param>
    public virtual void UpdateFractalDouble(double infx, double supx, double infy, double supy, int maxIt)
    {
        this.infx = infx;
        this.supx = supx;
        this.infy = infy;
        this.supy = supy;
        this.maxIteration = maxIt;

        latox = Math.Abs(supx - infx);
        latoy = Math.Abs(supy - infy);
        incx = latox / width;
        incy = latoy / height;

        this.Game.Window.Title = "Render WIP";
        complete = false;
        indici.Clear();
        for (int i = 0; i < total; i++)
        {
            int r = i / width;
            int c = i % width;
            indici.Enqueue(i);
            iterazioni[r, c] = 0;
            valori[r, c].set(0, 0);
            pixData[i] = Color.Lerp(pixData[i], Color.Black, 0.7f);
        }
        fractal.SetData(pixData);
    }

    public bool Complete
    {
        get { return complete; }
    }

    public int Width
    {
        get { return width; }
    }

    public int Heigth
    {
        get { return height; }
    }

    /// <summary>
    /// ridimenziona la texture su cui andiamo a calcolare il frattale
    /// </summary>
    /// <param name="width"></param>
    /// <param name="height"></param>
    public void resizeTexture(int width, int height)
    {
        //aggiorna i valori
        this.width = width;
        this.height = height;

        total = width * height;
        //ricrea il vettore dei colori
        pixData = new Color[total];

        //elimina la vecchia texture e la ricrea
        fractal.Dispose();
        fractal = new Texture2D(this.Game.GraphicsDevice,
            width,
            height, 1,
            TextureUsage.None, SurfaceFormat.Color);

        //ricrea la matrice con il numero di iterazioni fatte per ogni punto
        iterazioni = new int[height, width];
        //ricrea la matrice di numeri complessi per memorizzare il valore di ogni punto
        valori = new Cmplx[width, height];
        for (int i = 0; i < width; i++)
            for (int j = 0; j < height; j++)
                valori[i, j] = new Cmplx();

        //forzo l'aggiornamento
        UpdateFractalDouble(infx, supx, infy, supy, maxIteration);
    }

    #endregion
}

Preferisco soffermarmi sul metodo calcolaFrattale() e su quello che fa:
quello che facciamo è controllare quanti elementi sono presenti nella coda degli indici (punti) per cu idevo calcolare il prossimo passo della successione (in questo modo so quante estrazioni fare).
Quindi comincio con estrarre l'elemento in testa alla coda, e calcolati l'indice di riga e di colonna trovo il valore della costante C da usare nel calcolo della successione per quel punto.
Il valore di C cambia a seconda del punto che stiamo considerando.
Poi recuperiamo il valore precedente per Z, e controlliamo se in modulo è minore di 4, o se non abbiamo già eseguito per il quel punto il massimo numero di iterazioni.
Se il punto ha valore in modulo maggiore di 4 non fa parte dell'insieme, quindi va colorato in modo oppurtuno. Se invece ho fatto il massimo numero di iterazioni e il valore del punto è ancora minore di 4 lo considero come appartenente all'insieme di Madelbrot e lo coloro di nero.
Per la colorazione del punto, uso il rapporto tra il numero di iterazioni fatte prima di trovare che il punto diverge e il numero massimo di iterazioni come parametro nella funzione che trasforma un colore in scala HSV in RGB (dove il rapporto è usato come H, cioè Hue, cioè la tonalità del colore).
Cambiando il modo in cui si usa il numero di iterazioni fatte nel calcolo del colore del pixel possono ottenere scale di colore diverse.

Come potete vedere osservando il codice, ogni volta che eseguiamo il metodo calcolaFrattale, viene calcolato un passo per ogni punto che ancora non è stato considerato non appartenente all'insieme. Questo vuol dire che mano a mano che si procede nel calcolo avrò sempre meno punti su cui calcolare la successione, e quindi il calcolo tenderà a diventare più veloce!

NOTA: perchè per ogni esecuzione del metodo calcolo un solo passo? Non posso calcolare tutto il frattale in una volta sola?
Certo si poteva fare anche in quel modo, ma poi non avrei potuto renderizzare lo stato intermedio del calcolo: in parole semplici avrei dovuto aspettare che il calcolo fosse concluso per poter far vedere il frattale, mentre in questo modo lo posso renderizzare anche quando è in fase di costruzione. In questo modo l'attesa per la completazione del frattale sarà più sopportabile ^^

Diamo adesso una piccola spiegazione sull'implementazione del metodo UpdateFractal():
questo metodo non fa altro che aggiornare lo stato del componente con i valori passati, e fa ripartire il cacolo del frattale con i nuovi dati.
Adesso io volevo (ed ho poi fatto) un effetto di transizione tra un'immagine di un fratale e quella successiva, solo che non "sapevo" dove metterla e alla fine l'ho messa in questo metodo.
Analizziamo in dettaglio ciò:

           for (int i = 0; i < total; i++)
           {
               int r = i / width;
               int c = i % width;
               indici.Enqueue(i);
               iterazioni[r, c] = 0;
               valori[r, c].set(0, 0);
               pixData[i] = Color.Lerp(pixData[i], Color.Black, 0.7f);
           }

in questo pezzo di codice scorro il vettore che contiene il colore di ogni pixel dell'immagine. All'inizio impostavo tutti i pixel a nero, solo che così il passaggio tra un'immagine e la successiva era troppo netta (almeno per i miei gusti), e la cancellazione istantanea dell'immagine mi dava noia.
Quindi ho deciso di scurire tutta l'immagine facendo una semplice interpolazione tra il colore attuale de pixel ed il colore nero. Successivamente, quando il calcolo del frattale verrà riavviato, si andrà a lavorare di nuovo sulla stessa texture, colorando in modo diverso i pixel. Si avrà quindi un effetto di transizione tra un'immagine e l'altra.

Spero di essere stato chiaro fino ad ora, dato che adesso viene la parte più incasinata, e cioè gestire ed agire sul frattale tramite gli strumenti offerti da XNA.
Quello che voglio fare è riuscire, con il mouse, a selezionare una zona del frattale su cui fare zoom, quindi modificare gli estremi della zona del piano su cui stiamo calcolando la successione.
Voglio poi poter incrementare e drecrementare il numero massimo di iterazioni da fare come la risoluzione della texture su cui calcolo il frattale (più pixel ho più punti dovrò considerare, quindi avrò una maggiore precisione).

Invece di stare a schiaffarvi qui il codice dell'applicazione ve lo lascio la leggere direttamente dal sorgente, dato che la cosa che mi ero proposto di fare era spiegarvi come avevo implementanto il calcolo del frattale.
Però il codice sarà abbastanza commentanto, in questo modo non dovreste avere difficoltà!

Ecco il link al file
XNA-Mandelbrot.rar

Spero di essere stato il più chiaro possibile...
in caso non lo sia stato (cosa più che probabile XD) non esitate a chidere chiarimenti, a dirmi dove ho sbagliato, o anche a mandarmi a quel paese se volete!

Ciauz! Alla prossima da Odino!