Fraktale w C#

Nawiązując do artykułu Andrzeja Stasiewicza omawiającego tworzenie fraktala Mandelbrota (PCWK 7/2005) oraz do opisanej w poprzednim numerze implementacji liczb zespolonych w C#, przepiszemy algorytm tworzący ów piękny fraktal z użyciem utworzonej miesiąc temu struktury Complex. W ten sposób kod stanie się bardziej przejrzysty i ściślej związany z matematyczną definicją zbioru.

Nawiązując do artykułu Andrzeja Stasiewicza omawiającego tworzenie fraktala Mandelbrota (PCWK 7/2005) oraz do opisanej w poprzednim numerze implementacji liczb zespolonych w C#, przepiszemy algorytm tworzący ów piękny fraktal z użyciem utworzonej miesiąc temu struktury Complex. W ten sposób kod stanie się bardziej przejrzysty i ściślej związany z matematyczną definicją zbioru.

O rientację w tym, czym jest zbiór Mandelbrota, zdobyliśmy już trzy miesiące temu, czytając artykuł Andrzeja Stasiewicza. Teraz chciałbym przeformułować poznane wówczas definicje tak, żeby użyć liczb zespolonych, które zaimplementowaliśmy miesiąc temu (zob. artykuł "Liczby zespolone w C#" z wrześniowego numeru PCWK).

Definicja i algorytm

Zbiór Mandelbrota można rozumieć jako zbiór takich liczb zespolonych c, dla których każdy wyraz ciągu określonego w następujący sposób:

z0 = 0, z1 = z02 + c, zn = zn-12 + c

ma promień mniejszy od 2 (jak pamiętamy, promień liczby zespolonej to pierwiastek kwadratowy sumy kwadratów części rzeczywistej i urojonej). Innymi słowy, dla liczb c należących do tego zbioru spełniony jest warunek |zn| < 2 przy dowolnie dużej wartości indeksu n.

Rysunek 1. Cały zbiór Mandelbrota.

Rysunek 1. Cały zbiór Mandelbrota.

Jak wiadomo, liczby rzeczywiste mogą być umieszczone na jednej osi. Liczby zespolone wymagają płaszczyzny. Oś pionowa reprezentuje wówczas część urojoną, a pozioma - rzeczywistą. W ten sposób każdy punkt na płaszczyźnie odpowiada jednej liczbie zespolonej o wyznaczonej przez współrzędne punktu części rzeczywistej i urojonej. Chcąc graficznie zaprezentować zbiór Mandelbrota, możemy zaznaczyć czarnym kolorem te punkty, które odpowiadają liczbom zespolonym c, należącym do zbioru (rysunek 1). Pozostałe zostawiamy białe. Jeżeli chcemy narysować fraktal na powierzchni formy, to musimy przetestować tyle liczb zespolonych, ile pikseli zawiera obszar użytkownika tej formy.

Rysunek 2. Zbiór Mandelbrota wraz z zaznaczonymi kolorem punktami nienależącymi do zbioru dla zakresu liczb c od 0.2016+i0.541 do 0.2175+i0.5521.

Rysunek 2. Zbiór Mandelbrota wraz z zaznaczonymi kolorem punktami nienależącymi do zbioru dla zakresu liczb c od 0.2016+i0.541 do 0.2175+i0.5521.

Zanim zabierzemy się do pisania kodu, musimy rozwiązać jeden bardzo istotny problem. Zgodnie z definicją, żeby sprawdzić, czy liczba zespolona należy do zbioru Mandelbrota, trzeba wykonać nieskończoną liczbę iteracji, aby się upewnić, że każdy z wyrazów ciągu ma promień nie większy od 2. Tego nie potrafią nawet najszybsze komputery. Dlatego tworząc rysunek, ograniczymy liczbę iteracji do - powiedzmy - 255 i odwrócimy pytanie: będziemy szukać tych liczb c, o których wiemy, że nie należą do zbioru. Liczby c, dla których |z255| < 2 zaznaczamy na rysunku kolorem czarnym jako rokujące nadzieję, że są elementem zbioru Mandelbrota (czego oczywiście nie możemy być pewni, nie kontynuując testu), a pozostałe zaznaczamy na biało. Aby dodatkowo uatrakcyjnić rysunek, można zamiast nudnej bieli pokolorować te piksele, które nie należą do zbioru, uzależniając ich barwę od liczby iteracji, po których promień |zn| stał się większy od 2.

Nasz algorytm można zatem zapisać w następującej sekwencji instrukcji (por. przedstawiony na rysunku 3 schemat blokowy tego algorytmu):

1. Wybieramy wartość zmiennej c odpowiadającą kolejnemu pikselowi rysunku (część rzeczywista ustalana na podstawie współrzędnej poziomej, urojona - pionowej). Inicjujemy zmienną z wartością 0. Inicjujemy zmienną przechowującą liczbę iteracji indeks wartością 0.

2. Sprawdzamy, czy liczba iteracji przekroczyła 255. Jeżeli tak, rysujemy czarny piksel i przechodzimy do punktu (1).

3. Obliczamy nową wartość zmiennej z, przypisując jej z*z+c.

4. Sprawdzamy, czy Abs(z)>=2. Jeżeli warunek jest prawdziwy, kolorujemy piksel w zależności od liczby iteracji (wartość zmiennej indeks) i przechodzimy do punktu (1). Jeżeli warunek nie jest spełniony, zwiększamy wartość indeks o jeden iprzechodzimy do punktu (3).

5. Algorytm kończy się, gdy przetestowane zostaną wszystkie liczby c odpowiadające powierzchni rysunku.

Implementacja

Rysunek 3. Schemat blokowy odpowiadający przedstawionemu obok algorytmowi sprawdzania, czy liczba c odpowiadająca punktowi na rysunku należy do zbioru Mandelbrota.

Rysunek 3. Schemat blokowy odpowiadający przedstawionemu obok algorytmowi sprawdzania, czy liczba c odpowiadająca punktowi na rysunku należy do zbioru Mandelbrota.

Mając przemyślany algorytm, możemy przejść do jego implementacji. Tworząc projekt, wykorzystamy ponownie Visual C 2003. Nie powinno być jednak problemu z odtworzeniem poniższych czynności np. w C# Builderze lub nawet w Delphi 2005. Oczywiście to, że skorzystamy ze struktury Complex, nie zwiększy ani znacząco nie zmniejszy wydajności programu. Zwiększy się natomiast znacznie przejrzystość kodu i bardziej widoczny będzie jego związek z przytoczoną wyżej definicją. Można mieć nadzieję, że dzięki temu maleje groźba pomyłek logicznych.

Naciskamy [Ctrl Shift N], aby utworzyć nowy projekt. Zaznaczamy ikonę Windows Application z grupy Visual C# Projects. W pole Name wpisujemy nazwę projektu (proponuję "Mandelbrot") i klikamy OK. Już po chwili zobaczymy widok projektowania z podglądem formy. Interaktywność aplikacji ograniczymy do minimum. Mówiąc wprost, użytkownik nie będzie miał żadnego wpływu na rysowany fraktal. Będzie on obliczany dla ustalonego na stałe obszaru płaszczyzny zespolonej. Bez większego wysiłku jednak da się aplikację rozwinąć w taki sposób, żeby umożliwić przybliżanie wybranych myszą obszarów płaszczyzny (zob. kod aplikacji MandelbrotEksplorator umieszczonej na płycie dołączonej do czasopisma).

Rysunek 4. Zbiór nazwany na cześć Gastona Julii dla c=0.25 + i0.55

Rysunek 4. Zbiór nazwany na cześć Gastona Julii dla c=0.25 + i0.55

Do projektu włączamy plik Complex.cs, zawierający definicję typu Complex (zob. artykuł "Liczby zespolone w C#" w archiwum PCWK_PROG, umieszczonym na płycie dołączonej do czasopisma). W tym celu np. za pomocą Eksploratora Windows kopiujemy ów plik do katalogu projektu (prawdopodobnie Moje dokumenty/Visual Studio Projects/Mandelbrot). W oknie Solution Explorer z menu kontekstowego wybieramy podmenu Add, a z niego Add Existing Item, wskazujemy plik Complex.cs i klikamy Open. Zobaczymy Complex.cs wśród plików projektu. Musimy jeszcze uwzględnić przestrzeń nazw, w której zdefiniowana jest klasa Complex. Naciskamy klawisz [F7], aby przejść do edytora kodu, i na początku pliku dodajemy do zbioru wierszy rozpoczynających się od słowa kluczowego using jeszcze jeden, w którym deklarujemy użycie przestrzeni nazw Zespolone:

using Zespolone;

Wracamy do podglądu formy, naciskając kombinację [Shift F7]. W oknie Properties przechodzimy do listy zdarzeń (przycisk błyskawicy z podpowiedzią o treści Events), odnajdujemy zdarzenie Paint i klikamy dwukrotnie związane z nim pole edycyjne. Powstanie metoda zdarzeniowa Form1_Paint, w której umieszczamy poniższe polecenia:

private void Form1_Paint(object sender

System.Windows.Forms.PaintEventArgs e)

{

Text="Czekaj...";

//min,max - zakres badanej

//plaszczyzny zespolonej

Complex min=new Complex(-2.2,-1.2);

Complex max=new Complex(1.2,1.2);

Bitmap bufor=new Bitmap(

ClientSize.Width,ClientSize.Height);

//x,y - wspolrzedne pikseli

for (int x=0;x<bufor.Width;x++)

for (int y=0;y<bufor.Height;y++)

{

Complex c=min+new Complex(

x*(max.Real-min.Real)/bufor.Width

y*(max.Imag-min.Imag)/bufor.Height);

Complex z=0;

byte indeks=0;

for(;indeks<255;indeks++)

{

z=z*z+c;

if (z.Norm>=4) break;

}

if (indeks==255)

bufor.SetPixel(x,y,Color.Black);

else

bufor.SetPixel(x,y,Color.FromArgb(

255-indeks,255-indeks,255));

}

e.Graphics.DrawImage(bufor,0,0);

Text="Zbiór Mandelbrota";

}

W rezultacie wszystkie punkty odpowiadające liczbom zespolonym należącym do zbioru Mandelbrota (przynajmniej zgodnie z wiedzą, jaką dają testy o 255 iteracjach) zaznaczone zostaną kolorem czarnym. Punkty nienależące do zbioru kolorowane są w skali od białego do niebieskiego w zależności od szybkości, z jaką rozbiega się ciąg zn (rysunek 2).

Wyjaśnienia wymaga obecność zmiennej bufor. Dzięki niej nie rysujemy bezpośrednio na płótnie formy, co spowolniłoby działanie programu, a korzystamy z pamięci bufora poza pamięcią ekranu. Po przygotowaniu rysunku szybko kopiujemy odpowiedni fragment pamięci, wywołując metodę Graphics.DrawImage. Ta sztuczka znacznie skraca czas, po jakim zobaczymy gotowy obraz.

Niewielkim kosztem można zmodyfikować algorytm w taki sposób, żeby tworzyć zbiory Julii (rysunek 4). Wystarczy zmodyfikować inicjowanie zmiennych c i z. Jeżeli to zmienną z zwiążemy z przedstawionym na rysunku obszarem liczb zespolonych, natomiast c będzie dowolnym parametrem zespolonym, to uzyskamy zupełnie inny fraktal. Na jego postać w istotny sposób wpływa wartość stałej c - polecam przeprowadzenie kilku eksperymentów z dobieraniem różnych wartości tej liczby i rysowaniem odpowiadających im fraktali.

Aby uzyskać zbiór Julii, wystarczy w powyższym kodzie zmodyfikować jedynie dwa wiersze przypisujące wartości zmiennym c i z:

Complex c=new Complex(0.25,0.55);

Complex z=min+new Complex(

x*(max.Real-min.Real)/bufor.Width

y*(max.Imag-min.Imag)/bufor.Height);

Inne fraktale można tworzyć, zmieniając zależność rekurencyjną z punktu (3) naszego algorytmu. Można spróbować wyższych potęg z lub bardziej złożonych funkcji względem z i c. Należy jednak uwzględnić czas potrzebny na testowanie każdego punktu - użycie bardziej złożonej funkcji może go znacznie wydłużyć. Aby nie zmniejszać wydajności, można wówczas zredukować liczbę iteracji, ale to obniży z kolei jakość tworzonego rysunku..

<hr size=1 noshade>Klasyczne kolory Mandelbrota

Fraktal Mandelbrota (a właściwie punkty leżące blisko, ale nie należące już do zbioru) są często rysowane z użyciem charakterystycznej "ognistej" palety kolorów. Aby nasz program przykładowy (Mandelbrot lub MandelbrotEksplorator na płycie) mógł generować takie obrazy, trzeba wprowadzić w nim drobne modyfikacje. Najpierw musimy dodać do pól klasy tablicę int[,] Paleta przechowującą składowe kolorów w palecie. Następnie w konstruktorze głównej klasy aplikacji, czyli w metodzie Form1() dodajemy kod generujący wartości RGB dla każdego elementu tablicy. Poniższy kod to znacznie uproszczony algorytm konwersji przestrzeni kolorów HSV na RGB dostosowany do generowania kolejnych barw spektrum światła widzialnego od czerwieni do fioletu.

Paleta = new int[256,3];

int[] GRB = new int[3] {0,255,0};

int d = 1;

for (int i=0; i<4; i++) {

if (GRB[i % 3] == 255) d=-1; else d=1;

for (int j=0; j<50; j++){

GRB[i % 3] = GRB[i % 3] + d*5;

Paleta[i*51+j,0] = (byte)GRB[1];

Paleta[i*51+j,1] = (byte)GRB[0];

Paleta[i*51+j,2] = (byte)GRB[2];

}

}

Zmieniona kolejność składowych w zmiennej GRB (zamiast RGB) umożliwia wykorzystanie funkcji modulo do wybierania kolejnej składowej, która będzie modyfikowania (stopniowo zwiększana lub zmniejszana) w zależności od parametru d, oznaczającego kierunek zmian. Na koniec przy wpisywaniu wartości do palety ponownie zmieniamy kolejność, aby powrócić do układu RGB. Na koniec trzeba zmodyfikować kod rysujący punkty, tak by korzystał z palety kolorów. Zmieniamy zatem linię wywołującą metodę SetPixel w procedurze obsługi zdarzenia Form1_Paint i gotowe - efekt widoczny jest na rysunku.

bufor.SetPixel(x,y,Color.FromArgb(Paleta[

indeks,0],Paleta[indeks,1],Paleta[indeks,2]);

Struktura Complex

W zeszłym miesiącu w artykule "Liczby zespolone w C#" pojawił się błąd w jednej z metod struktury Complex. Funkcja Sqrt powinna wyliczać pierwiastek, korzystając z modułu liczby zespolonej, a więc z.Radius a nie z.Norm. Oto prawidłowa postać funkcji:

public static Complex Sqrt(Complex z) {

return Math.Sqrt(z.Radius)*(new Complex(Math.Cos(

0.5*z.Phase),Math.Sin(0.5*z.Phase)));

}


Zobacz również