Fraktal Mandelbrota

Dziś skonstruujemy zbiór Mandelbrota a także trochę pokolorujemy (odcieniami szarości) jego otoczenie. Postaram się wyjaśnić to w szczegółach. Zawodowy programista będzie się z pewnoscią znów zżymał, że albo zbyt szczegółowo, albo niedostatecznie szczegółowo.

Ostatnim razem zajmowaliśmy się iteracjami transformacji

T(z) = z²+c

gdzie z i c są liczbami zespolonymi. Startowaliśmy od z=0 i wykonywaliśmy tą transformację wiele razy, przy ustalonym c. Interesowało nas przy tym za którym razem nastapi ucieczka z „zaczarowanego koła” o promieniu 2, tzn dla jakiego n, T(T(T(....(T(z))))) (n razy) otrzymana w ten sposób liczba zespolona będzie miała kwadrat modułu większy od 4.

Interesowaliśmy się przy tym wyłącznie stałymi c postaci:

c = -0.75+0.1j

c = -0.75+0.01j

...

c = -075+0.0000001j

Zauważyliśmy ciekawą prawidłowośc – pojawiła się ni stąd ni z owąd liczba Pi.

Dziś zbadamy większy zakres stałych c. Część rzeczywistą c będziemy testować od -2.5 do 1.5,

zaś część urojoną będziemy testować w zakresie od -2 do 2.

Czemu akurat zakres zmienności stałej c wybierzemy? Po prostu dlatego, że, jak się okaże, jest on interesujący. Tak czy siak interesuje nas jednynie wnętrze zaczarowanego koła. Mógłbym się więc ograniczyć do części rzeczywistej zaczynającej się od -2, po co nam -2.5, które od początku jest poza kołem? Powód jest prosty: chcę mieć kwadrat o długości 4 i wysokości 4 i chcę mieć obrazek fraktala w miarę wycentrowany w środku tego kwadratu. Kwestia prób i błędów. Każdy będzie mógł to sobie zmienić, kiedy metoda rysowania fraktala stanie się jasna. Za chwilę to wyjaśnię, muszę jednak zacząć od wyjaśnienia jak będziemy programować grafikę. Ogranicze się przy tym do grafiki z FreeBasic.

Okno graficzne składa się z pikseli. Ileś tam pikseli w poziomie i ileś tam pikseli w pionie. Powiedzmy pw pikseli w poziomie i ph pikseli w pionie (pw od pixel-width, ph od pixel – height). W programie, którym będziemy się bawić wezmę pw=440 i ph=440. Dlaczego 440? Bo to jest bezpieczna szerokośc grafiki w salonie24. Czytelnik będzie mógł to u siebie zmienić bacząc jednak, by nie przekroczyć rozdzielczości ekranu – inaczej komputer może się zezłościć i nastąpi krach systemu – reset. Trzeba przy tym pamiętać, że czym większa rozdzielczośc, tym dłużej będzie trwało wykonywanie programu, bo tym więcej pikseli trzeba będzie obsłużyć, a dla obsłużenia każdego piksela potrzebne będą rachunki. Dla jednych krótkie, dla innych długie.

Najpierw kompletny tekst programu, tym razem bez komentarzy wewnątrz, po czym wyjaśnienie jego działania.

''''''

#define stopat 500
#define pw 440
#define ph 440

screenres pw,ph,32

Declare Sub drawfractal


drawfractal
Sub drawfractal
Dim As Integer i, i_out, gray
Dim As Integer m,n
Dim As Double cr,ci,x,y,pom,a
Dim As Double delta_h=4/pw,delta_v=-4/ph
Dim As Double c0_x=-2.5+2/pw, c0_y=2-2/pw
Dim As Integer Ptr sp=screenptr
screenlock
For n=0 To ph-1
For m=0 To pw-1
cr=m*delta_h+c0_x
ci=n*delta_v+c0_y
x=0
y=0
i=0
i_out=stopat
Do While i< stopat
i=i+1
pom = x*x-y*y+cr
y=2*x*y+ci
x=pom
a=x*x+y*y
if a>4 then
i_out=i
Exit Do
end if
Loop
If i_out
gray=Abs((i_out And 63)-31)*7+31
*(sp+n*pw+m)=rgb(gray,gray,gray)
Else
' *(sp+n*pw+m)=0
End If
Next
Next
screenunlock
Sleep
End Sub

W programie wystąpią dwie zmienne przyjmujace wartości całkowite m,n. Zmienna m będzie się zmieniać w przedziale od 0 do pw-1, zmienna n od 0 do ph-1.

m=0, n=0 to piksel w lewym górnym rogu okna graficznego

m=pw-1, n=ph-1 to piksel w prawym dolnym rogu okna

Będziemy szli wierszami od góry do dołu:

n=0, m=0,1,2,...,pw-1 to pierwszy wiersz od góry

n=1, m=0,1,2,..., pw-1 to drugi wiersz od góry

i tak dalej, aż do

n=ph-1, m=0,1,2,..., pw-1 to ostatni, dolny wiersz okna.

Teraz musimy przetłumaczyć każdy piksel na liczbę zespoloną c reprezentowaną przez ten piksel. Jak to robimy? Potrzebny jest prosty rachunek. Mamy pw pikseli w poziomie a szerokośc okna w liczbach wybraliśmy równą 4. Zatem na każdy piksel przypada delta_x = 4/pw osi x. Jako, że wybraliśmy kwadrat 4x4, będzie na każdy piksel przypadało delta_y = 4/ph osi y.

Przystąpię teraz do opisu samego programu rysującego fraktal. Zaczyna się od deklaracji stałych. Poprzednio deklarowaliśmy przez „const”. Można jednak je zadeklarować i tak:

#define stopat 500

#define pw 440

#define ph 440

Stała stopat deklaruje maksymalna liczbę iteracji transformacji T. Mówi ona o tym, ile razy będziemy powtarzać transformację T by się przekonać czy punkt nie ucieknie z koła. Punkty, które nie uciekną do 500 razy będziemy uważać za „nie uciekną nigdy” i będziemy malować na czarno. Te będą tworzyć nasz fraktal. Troche to ryzykowne, nieprawdaż? Ryzykowne, ale działa. Warto w domu wypróbować najpierw stopat 5, potem stopat 50, wreszcie stopat 500. Można i wyjść poza 500, wtedy czas się odpowiednio wydłuży. Ale nie tylko to, bo zacznie odgrywac role skończona precyzja obliczeń komputera. Na obliczeniach w których podnosimy liczby do potęgi 1000 nie można z całym sercem polegać.

W następnej linii programu deklarujemy rozdzielczość ekranu i głębię koloru:

screenres pw,ph,32

Co to jest pw i ph już wiemy. Liczba 32 oznacza, że będziemy uzywać 32-bitowego TrueColor. W praktyce oznacza to, że każdy z kolorów Red, Green, Blue będzie mógł przyjmować wartość od 0 do 253. Dzisiaj chyba wszystkie karty graficzne taką głębią operują. Kolorować (zresztą jedynie odcieniami szarości) będziemy dopiero na końcu.

Jedna uwaga: poniewaź program będzie się dobierał do pamięci komputera, istnieje zawsze niebezpieczeństwo, że coś będzie nie tak. Czyjś komputer będzie miał np. mniejszą głębię kolorów niż zakładana. Trudno przwidzieć co się w takich przypadkach stanie. Dlatego przed uruchomieniem po raz pierwszy programu trzeba koniecznie zapisać na dysk wszystko wszystko czego nie chcielibyśmy stracić w przypadku, gdyby komputer nagle, w czasie wykonywania programu, się zresetował. Lepiej być przygotowanym na najgorsze i cieszyć się, gdy nic złego się nie stanie niż na odwrót: najpierw się cieszyć, a potem płakać, bo straciliśmy niedokończoną inna pracę. Prawdę mówiąc dotyczy to każdego innego uruchomianego po raz pierwszy programu, również tych pochodzących z Mikrosoftu!

Samo obliczanie i kolorowanie zawrzemy w podprogramie. Nie jest to konieczne tym razem, jednak przyda się w przyszłości, kiedy nauczymy się robić ZOOM na wybrane otoczenie przy pomocy kliknięcia myszką. Zatem z myślą o przyszłości deklarujemy podprogram, z angielska - subrytynę

Declare Sub drawfractal

Nazałem ją też z angielska – drawfractal. Tutaj nie potrzebuje ona żadnych parametrów. Przy zoomie będzie miała parametry punktu w którym klikniemy myszą.

Dalej wywołujemy podprogram, każemy mu działać:

drawfractal

A następnie opisujemy na czym to działanie polega:

Sub drawfractal

.....

End Sub

Zajmijmy się teraz tym co jest w środku, pomiędzy „Sub drawfractal” a „End Sub”. Tam właśnie odbywa się właściwa praca. Zaczynamy znów od deklaracji zmiennych występujacych w podprogramie:

Dim As Integer i, i_out, gray

Dim As Integer m,n

Dim As Double cr,ci,x,y,pom,a

Dim As Double delta_h=4/pw,c0_x=-2.5+2/pw, delta_v=-4/ph, c0_y=2-2/pw

Dim As Integer Ptr sp=screenptr

Zmienna i będzie przyjmować wartości całkowite i numerować iteracje.

Zmienna i_out, również całkowita, będzie notować przy której iteracji punkt ucieka z koła.

Zmienna gray, całkowita, będzie przypisywała odcień szarości pikselowi, zależnie od wartości i_out.

Czym szyciej ucieka, tym jasniejszy piksel. Nie ucieknie aż do stopat – piksel czarny. Ale to dopiero pod koniec podprogramu.

Zmienne m,n – całkowite będą numerować wiersze i kolumny okna – jak to wyjasniałem wyżej.

Zmienne cr, ci będą liczbami rzeczywistymi w podwójnej precyzji (16 miejsc znaczących) i będą służyły do przechowywania części rzeczywistej i urojonej stałej c. Zmienne x,y będą służyły do przechowywania i uaktualniania rzeczywistej i urojonej współrzędnej transformowanego punktu z.

Zmienna pom będzie służyć do chwilowego przechowania rzeczywistej wartości x, zanim nie policzymy części urojonej. Tak jak w poprzedniej nocie.

Dalej mamy linię z

Dim As Double delta_h = 4/pw, cr = -2.5+2/pw, delta_v = -4/ph, ci = 2-2/pw

Tutaj deklarujemy dalsze zmienne jednocześnie nadając im wartości.

delta_h=4/pw oraz delta_v=-4/ph już omówiłem. Dlaczego jednak pojawił się minus w wyrażeniu na delta_v? Idzie o to, że okno graficzne w pionie numeruje się od góry do dołu. Zaś współrzędna y, normalnie, zwiększa się od dołu do góry. Zwiększając zatem numer wiersza zmniejszamy wartośc y-ka. Stąd i minus.

Pozostają do wyjasnienia tajemnicze: c0_x=-2.5+2/pw oraz c0_y=2-2/pw. Są to kartezjańskie wspólrzędne przypisywane lewemu górnemu pikselowi. Lewy górny róg ekranu ma mieć wspólrzędne

(-2,5,2). To lewy górny róg piksela. Prawy dolny róg tego piksela będzie miał współrzędne (-2.5+delta_x, 2-delta_y). No to proszę policzyć, używając wartości delta_w=4/pw i delta_h=4/ph, jakie wspólrzędne będzie miał środek tego piksela? Ta niewielka poprawka powinna poprawiać precyzję oddania kształtu fraktala w skali jednego piksela. Trzeba pamiętać o tym, że pewne punkty są szczególnie czułe na zmianę wartości stałej c. Do takich należy punkt padany w poprzedniej nocie, w otoczeniu c=-075+0j. Widzieliśmy, że zmana części urojonej o jedną setną zwiększa ilość iteracji potrzebnych do ucieczki stukrotnie. Punkty szczególnie czułe, to punkty na brzegu fraktala, brzeg ten zaś jest postrzępiony – fraktalny!

Pozostaje do wyjasnienia linia:

Dim As Integer Ptr sp=screenptr

Nie chcę tu wchodzić w szczegóły. Ta linia po prostu rezerwuje w pamięci RAM komputera obszar potrzebny do przechowania zawartości ekranu zanim ten obszar, po „narysowaniu go w pamięci” zostanie jedną komendą wrzucony cały na ekran. Zamiast tego moglibyśmy rysować na prawdziwym ekranie, piksel za pikselem. Spowolniałowby to jednak działanie programu. Ta linia programu jest taka sama cokolwiek byśmy nie rachowali. Zobaczymy w dalszym ciągu jak takie rysowanie w pamięci się odbywa.

Najpierw zamrażamy okno graficzne komendą:

screenlock

po czym przystepujemy do obróbki pikseli, jeden po drugim, wierszami od lewo-góra do prawo dół.

Służy do tego podwójna pętla:

For n=0 To ph-1

For m=0 To pw-1

........

Next

Next

Po zakończeniu obróbki wszystkich pikseli okna odmrażamy ekran, czyli wrzucamy z pamięci RAM do pamięci graficznej – obraz ukazuje się na ekranie komputera po wykonaniu komendy

screenunlock

Po czym zamrażamy ten obraz, komendą

Sleep

by natychmiast po ukazaniu się nie zniknął. Takie postępowanie jest specyficzne dla tego dialektu Basica. W Pascalu będzie to wyglądało inaczej, zresztą też zależy od wersji, od szczególnej implementacji danego języka programowania.

Zajmijmy się teraz tym, co jest we wnętrzu pętli. Przy danych n,m obrabiamy piksel w m-tej kolumnie i n-tym wierszu ekranu. Jakie są współrzędne x,y środka tego piksela. Każdy może się przekonać, robiąc prosty szkic ekranu, numeracji pikseli i współrzędnych brzegu ekranu, że współrzędne te wynoszą:

cr=m*delta_h+c0_x

ci=n*delta_v+c0_y

I to jest wartość stałej c przypiywana pikselowi (m,n). Dalej już idziemy, z niewielkimi tylko zmianami tak jak w poprzednim programie. Nadajemy najpierw naszym zmiennym x,y,i wartości poczatkowe:

x=0

y=0

i=0

Pierwsze dwie linijki oznaczają, że naszym iteracjom bedzie za każdym razem podlegał punkt z=0, zaś i=0 oznacza, że na poczatku mamy zero iteracji. Mamy jeszcze jedną linię inicjującą:

i_out = stopat.

Nadajemy poczatkowo zmiennej i_out największą dopuszczalną wartość. Oznacza to, że na samym początku traktujemy kazdy piksel jako czarny – należący do fraktala. Następnie rozpoczynamy pętle iteracji i jesli punkt ucieknie wcześniej, wtedy zmienimy kolor obrabianego piksela. Czym szybciej ucieknie, tym jaśniejszy kolor.

Dalej zaczyna się petla, prawie tak samo jak w poprzednim programie, z małymi jednak zmianami:

Do While i< stopat

i=i+1

pom = x*x-y*y+cr

y=2*x*y+ci

x=pom

a=x*x+y*y

if a>4 then

i_out=i

Exit Do

end if

Loop

Przede wszystkim petla zaczyna się od

Do While i< stopat

i=i+1

miast, jak poprzednio, tylko

Do

i=i+1

Dlaczego? Otóż w poprzednim programie popełniłem wykroczenie za które dyplomowany programista powinien być karany, przynajmniej upomnieniem ustnym. Każemy bowie i zwiększać się o jeden, nie mając pewności, że temu zwiększaniu się nastapi kiedyś kres. Jeśli taki kres nie nastąpi, to w pewnym momencie wartośc i przekroczy największą dopuszczlną dla danego komputera i typu zmiennej liczbę całkowitą i nasz program się „wywali”. Albo też, jeśli obliczanie jest wolne, po godzinie pracy będziemy musieli zresetować komputer nie otrzymawszy końcowego wyniku. Tym razem zabezpieczamy się z góry. Jak powiedziałem wyżej, 500 iteracji będziemy traktować jako „praktyczną nieskończoność”. Nic się nie zdarzyło do 500? Piksel czarny. Punkt nie uciekł. Wartośc stałej c należy do fraktala. Następny piksel. I tak dalej.

Dalej leci jak w poprzednim programie, z mała tylko zmianą:

Do While i< stopat

i=i+1

pom = x*x-y*y+cr

y=2*x*y+ci

x=pom

a=x*x+y*y

if a>4 then

i_out=i

Exit Do

end if

Loop

Gdy tylko wyskakujemy poza koło, notujemy przy której iteracji wyskakujemy – przechowujemy tą liczbę w i_out, przerywamy dalsze iterowanie. Pozostaje zabarwienie odpowiedniego piksela zależnie od wartości i_out. Temu służy fragment:

If i_out<>stopat Then

gray=Abs((i_out And 63)-31)*7+31

*(sp+n*pw+m)=rgb(gray,gray,gray)

Else

*(sp+n*pw+m)=0

End If

Ten fragment nie jest istotny dla poczatkujących. Ci mogą poniższe wyjasnienie tylko przelecieć. Prawdę mówiąć, z ręką na sercu, nie jestem w 100% pewien czy sam go rozumiem. Specjaliści mnie poprawią, jeśli coś nakłamię. Ideę tego fragmentu zaczerpnąłem z innego programu, spróbowałem zrozumieć i zdecydowałem się zostawić jak jest, bez zmian. Kto zechce, kto potrafi, ten tutaj poeksperymentuje. Zapewne wymyśli lepszą alternatywę. Zacznę od tego, że nie jest to tak, że czym większe i_out, tym ciemniejszy piksel. Nie byłoby to wskazane. Chcemy bowiem dość ostro zobaczyć granice fraktala. Przy granicy i_out zmienia się szybko. Zaś oko ludzkie nie odróżni całkiem czarnego od prawie czarnego. Granice frakatala zatem rozmyłyby się. Zatem zamiast stopniowego ściemniania, w miarę zbliżania się do granicy fraktala, stosujemy tu cyrkulację 33 odcieni szarości. Są one dane liczbami:

31,38,45,52,59,66,73,80,87,94,101,108,115,122,129,136,143,150,157,164,171,178,185,192,199,206,213,220,227,234,241,248,255

Jak widać zmiany następują o siedem. 255 to kolor biały. 0 -zarezrewowane dla punktów które nie uciekają aż do stopat, to kolor cazrny. Ze wzrostem i_out o jeden przechodzimy do następnego odcienia szarości z listy. Kiedy listę wyczerpiemy, wracamy do jej poczatku. Załatwia to formuła

gray=Abs((i_out And 63)-31)*7+31

Nie będę w nią wchodził, zwrócę tu tylko uwagę, że występuje w niej „And” operujące na bitach i zastępujace tu funkcję „minimum”. Wązne jest w tej operacji to, że 63 w zapisie dwójkowym to 111111 – same jedynki. Tyle o przypisywaniu odcieni szarości pikselom zależnie od wartości i_out.

Wreszcie ostatnia linia:

*(sp+n*pw+m)=rgb(gray,gray,gray)

Tutaj sp to adres początku pamięci zarezerwowanej dla naszego obrazka w pamięci komputera. Nie wchodząc w szczegóły, sp+n*pw+m to adres piksela w m-tej kolumnie i w n-tym wierszu. A cała linijka zapisuje pod ten adres kolor szary – tyle samo czerwonego co zielonego co niebieskiego. Specjalista z pewnością to uściśli, jesli poczuje taka potrzebę.

Tyle na dziś. Zyczę powodzenia, uruchomienia programu, zobaczenia czegoś jak na obrazku poniżej. Potem można poeksperymentować zmieniając stopat i zmieniając pw i ph – byle nie przekroczyć rozdzielczości ekranu na naszym komputerze.

Ach, jeszcze jedno: sam program trzeba skopiować do edytora FreeBasic i zapisać na dysku nadając mu nazwę. Ja nadałem swojemu nazwę „mymandel2.bas”:

Czy da się coś podobnego zrobić w Excelu? Nie jestem pewien ....