Buddhabrot

[Ashish]

Implementation of https://en.wikipedia.org/wiki/Buddhabrot in Qb64.

Here's the code -

Code: QB64: [Select]

1. 'Buddhabrot Fractal By Ashish 12 August, 2018
2. 'https://en.wikipedia.org/wiki/Buddhabrot
3.  
4.  
5. _TITLE "Buddhabrot Fractal"
6.  
7. SCREEN _NEWIMAGE(600, 600, 32)
8. DIM r_plot(_WIDTH, _HEIGHT), g_plot(_WIDTH, _HEIGHT), b_plot(_WIDTH, _HEIGHT) 'for red,green and blue channels
9. DIM SHARED seq_x(5000), seq_y(5000)
10.  
11. RANDOMIZE TIMER
12.  
13. PRINT "Generating density plot"
14. PRINT
15. max_itr = 1000000
16. t# = TIMER
17.  
18. 'for speed
19. $CHECKING:OFF
20. WHILE itr < max_itr
21.     FOR b = 0 TO 10
22.         x = RND * _WIDTH
23.         y = RND * _HEIGHT
24.  
25.         xx = map(x, 0, _WIDTH, -2, 2)
26.         yy = map(y, 0, _HEIGHT, -2, 2)
27.  
28.         n = num_itr(xx, yy, 170) 'only those points, which escape to infinity are useful.
29.  
30.         'we are now going to store density of each series of this point.
31.         'I've used 3 arrays for color and in each array, the density is stored for different number of iterations.
32.         IF n > 1 THEN
33.             FOR j = 0 TO n - 2
34.                 px = ABS(INT(map(seq_x(j), -2, 2, 1, _WIDTH - 1)))
35.                 py = ABS(INT(map(seq_y(j), -2, 2, 1, _HEIGHT - 1)))
36.                 IF px > 0 AND px < _WIDTH AND py > 0 AND py < _HEIGHT THEN
37.                     r_plot(px, py) = r_plot(px, py) + 1
38.                     IF n < 40 THEN b_plot(px, py) = b_plot(px, py) + 1
39.                     IF n < 100 THEN g_plot(px, py) = g_plot(px, py) + 1
40.                 END IF
41.             NEXT
42.         END IF
43.     NEXT
44.     itr = itr + 1
45.     IF itr MOD 10000 = 0 THEN PRINT ".";
46. WEND
47. tt# = TIMER - t#
48. $CHECKING:ON
49. 'NEXT x, y
50.  
51. CLS
52. FOR y = 0 TO _HEIGHT
53.     FOR x = 0 TO _WIDTH
54.         PSET (y, x), _RGB(r_plot(x, y), g_plot(x, y), b_plot(x, y))
55. NEXT x, y
56. PRINT tt#; "s"
57.  
58.  
59. FUNCTION num_itr (x_0, y_0, max_itr)
60.     x = 0: y = 0
61.     x_n_1 = 0: y_n_1 = 0
62.  
63.     FOR i = 0 TO max_itr
64.         x_n_1 = x * x - y * y + x_0
65.         y_n_1 = 2 * x * y + y_0
66.         x = x_n_1
67.         y = y_n_1
68.         seq_x(i) = x: seq_y(i) = y
69.         IF ABS(x * x + y * y) > 4 THEN num_itr = i: EXIT FUNCTION
70.     NEXT
71. END FUNCTION
72.  
73. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
74.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
75. END FUNCTION
76.  
77.  

 This image is generated with this code of about 50000000000 cycles.

This image is also generated with this code of about 1700000000 cycles

[bplus]

Hi Ashish,

That is a very nice version you have made!

You might enjoy comparing to this code I found somewhere:

Code: QB64: [Select]

1. 'Buddhabrot by Antoni Gual 6/2003 agual@eic.ictnet.esB
2. '------------------------------------------------------------------------------
3. 'This different rendering of the Mandelbrot Set, created by Melinda Green
4. ' reminds the shape of the classic sculpture of the sitting Buddha.
5. '-----------------------------------------------------------------------------
6. 'The number of calculations needed to obtain this shape is much bigger than
7. ' for classic mandelbrot set.The complete drawing can take hours
8. 'Good news is you can divide time needed by 5
9. ' if you load qb with ffix library, a patch by V1ctor for the QB slow floting
10. ' point operations. You can find ffix in my site among many others.
11. '    www.geocities.com/antonigual/qbasic
12. '
13. ' You must uncomment the two lines invokking ffix,then call qb with:
14. '  c:\qbpath\qb budabrot.bas /lc:\qbpath\ffix
15. '-----------------------------------------------------------------------------
16. 'Classic Mandelbrot is drawn by plotting the points in the complex plane where
17. 'the Mandelbrot formula leads to a convergent series.
18. '
19. 'The BuddhaBrot is drawing by cumulating all intermediate points of the
20. 'diverting Mandelbrot series.
21. '----------------------------------------------------------------------------
22. 'Two sites about Mandelbrot and its variants:
23. 'http://mandelbrot.dazibao.free.fr/                     'Mandelbrot Set in QB
24. 'http://www.superliminal.com/fractals/bbrot/bbrot.htm   'Melinda Green's site
25. '
26. 'Have Fun!...E-mail me yor improvements on this program..
27. '------------------------------------------------------------------------------
28. 'If you are using ffix, uncomment the next two lines to boost speed x 5
29. ''DECLARE SUB ffix
30. ''ffix
31.  
32. SCREEN 13
33.  
34. 'set palette made of shades of green
35. OUT &H3C8, 0
36. FOR c = 0 TO 254
37.     temp = c / 4
38.     OUT &H3C9, 0
39.     OUT &H3C9, temp
40.     OUT &H3C9, 0
41. NEXT
42.  
43. 'plus awhite color for the time display
44. OUT &H3C9, 63
45. OUT &H3C9, 63
46. OUT &H3C9, 63
47. COLOR 255
48.  
49. CLS
50.  
51. 'If you increase niter you will have a crisper but slower rendering
52. CONST niter = 300
53.  
54. 'set up arrays to save intermediate points
55. DIM real(niter)
56. DIM imag(niter)
57. t! = TIMER
58.  
59. DO
60.     'seed the mandelbrot series with random points of complex plane
61.     x = 4 * RND - 2
62.     y = 4 * RND - 2
63.     im = 0: re = 0: re2 = 0: im2 = 0
64.     diverts% = 0
65.     'iterate the mandelbrot formula
66.     FOR iter% = 0 TO niter
67.         im = 2 * re * im + x
68.         re = re2 - im2 + y
69.         'Save the intermediate points of the series to arrays
70.         'IF re > 1E+21 THEN EXIT FOR
71.         'IF im > 1E+21 THEN EXIT FOR
72.         real(iter%) = re
73.         imag(iter%) = im
74.         im2 = im * im
75.         re2 = re * re
76.         'If the seed point leads to a diverting series, we are interested
77.         IF re2 + im2 > 4 THEN diverts% = 1: EXIT FOR
78.     NEXT
79.  
80.     'Plot to screen the points of the DIVERGENT series
81.     IF diverts% THEN
82.         FOR it% = 0 TO iter%
83.             xs% = CINT(imag(it%) * 80) + 160
84.             ys% = CINT(real(it%) * 80) + 120
85.  
86.             pp% = POINT(xs%, ys%)
87.             IF pp% < 255 THEN PSET (xs%, ys%), pp% + 1
88.         NEXT
89.         'To display rendering time, not really needed
90.         np% = (np% + 1) AND &HFF
91.         IF np% = 0 THEN LOCATE , 1: PRINT USING "####.#"; TIMER - t!;
92.     END IF
93. LOOP UNTIL LEN(INKEY$)
94.  

[_vince]

neat

[FellippeHeitor]

Math is a thing of beauty, right?

[Ashish]

Thank you Bplus, Vincent and Fellippe!
@bplus
That code is quite good. I like the way it renders in real time. Following this, I edited my code so that it can also
render the fractal in real time. Now, you are also allowed to enter you own number of iteration per pixel.

Code: QB64: [Select]

1. 'Buddhabrot Fractal By Ashish 12 August, 2018
2. 'https://en.wikipedia.org/wiki/Buddhabrot
3. 'EDITED : 13 August, 2018
4. 'Now renders in real time, as well as allow you to choose number of iterations per pixel.
5.  
6.  
7. _TITLE "Buddhabrot Fractal"
8.  
9. SCREEN _NEWIMAGE(600,600, 32)
10. DIM r_plot(_WIDTH, _HEIGHT), g_plot(_WIDTH, _HEIGHT), b_plot(_WIDTH, _HEIGHT) 'for red,green and blue channels
11.  
12. RANDOMIZE TIMER
13.  
14. PRINT "Enter the maximum iteration for each pixel : "
15. INPUT temp
16. DIM SHARED seq_x(temp), seq_y(temp)
17.  
18. rf! = RND: gf! = RND: bf! = RND
19.  
20.  
21. max_itr = 1000000
22. t# = TIMER
23.  
24. 'for speed
25. $CHECKING:OFF
26. WHILE itr < max_itr
27.     FOR b = 0 TO 15
28.         x = RND * _WIDTH
29.         y = RND * _HEIGHT
30.  
31.         xx = map(x, 0, _WIDTH, -2, 2)
32.         yy = map(y, 0, _HEIGHT, -2, 2)
33.  
34.         n = num_itr(xx, yy, temp) 'only those points, which escape to infinity are useful.
35.  
36.         'we are now going to store density of each series of this point.
37.         'I've used 3 arrays for color and in each array, the density is stored for different number of iterations.
38.         IF n > 1 THEN
39.             FOR j = 0 TO n - 2
40.                 px = ABS(INT(map(seq_x(j), -2, 2, 1, _WIDTH - 1)))
41.                 py = ABS(INT(map(seq_y(j), -2, 2, 1, _HEIGHT - 1)))
42.                 IF px > 0 AND px < _WIDTH AND py > 0 AND py < _HEIGHT THEN
43.                     IF n < temp * rf! THEN r_plot(px, py) = r_plot(px, py) + 1
44.                     IF n < temp * bf! THEN b_plot(px, py) = b_plot(px, py) + 1
45.                     IF n < temp * gf! THEN g_plot(px, py) = g_plot(px, py) + 1
46.                     PSET (py, px), _RGB(0, 0, 0)
47.                     PSET (py, px), _RGB(r_plot(px, py), g_plot(px, py), b_plot(px, py))
48.                 END IF
49.             NEXT
50.         END IF
51.     NEXT
52.     itr = itr + 1
53.     IF itr MOD 10000 = 0 THEN _DISPLAY
54. WEND
55.  
56.  
57. FUNCTION num_itr (x_0, y_0, max_itr)
58.     x = 0: y = 0
59.     x_n_1 = 0: y_n_1 = 0
60.  
61.     FOR i = 0 TO max_itr
62.         x_n_1 = x * x - y * y + x_0
63.         y_n_1 = 2 * x * y + y_0
64.         x = x_n_1
65.         y = y_n_1
66.         seq_x(i) = x: seq_y(i) = y
67.         IF ABS(x * x + y * y) > 4 THEN num_itr = i: EXIT FUNCTION
68.     NEXT
69. END FUNCTION
70.  
71. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
72.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
73. END FUNCTION
74.  
75.  
76.  

Image is generated with 20 iterations per pixel. As the number of iterations increases, the fractal
becomes more smooth and diverse.

[bplus]

Nice mod Ashish,

It sure beats watching dots go across screen for awhile.

I recommend CLS after input, leave checking alone specially now that it draws in real time, you don't need the speed so much because it's nice to watch image develop.

I like how there is a color change from just one increment if you pick a low iteration number, eg 20, 25, 26 so each run is unique.

There is ideal range for numbers eg 1000 runs too long after a certain whiteness is reached and 10000 just seems too slow  without anything really new or significant developing (which was when I was looking to check out but Alt+F4 did not work!).

All in all, I'd mark this you "best answer" so far.

[_vince]

Yeah, nice mod Ashish.  Do you think you can make it so that you could zoom in and see the smaller details in this beautiful fractal?

[Ashish]

There is ideal range for numbers eg 1000 runs too long after a certain whiteness is reached and 10000 just seems too slow  without anything really new or significant developing (which was when I was looking to check out but Alt+F4 did not work!).

Buddhabrot can take hour to render. ;)
See the image on wikipedia having 20,000 and 1,000,000 iterations per pixel.


@Vincent
Now, the program support zooming feature. Zoom into the wonderful world of Buddhabrot!

Code: QB64: [Select]

1. 'Buddhabrot Fractal By Ashish 12 August, 2018
2. 'https://en.wikipedia.org/wiki/Buddhabrot
3. 'EDITED : 13 August, 2018
4. 'Now renders in real time, as well as allow you to choose number of iterations per pixel.
5. 'EDITED : 14 Augist, 2018
6. 'Added Zooming Feature
7.  
8.  
9. _TITLE "Buddhabrot Fractal"
10.  
11. TYPE vec2
12.     x AS SINGLE
13.     y AS SINGLE
14. END TYPE
15.  
16. SCREEN _NEWIMAGE(500, 500, 32)
17. DIM r_plot(_WIDTH, _HEIGHT), g_plot(_WIDTH, _HEIGHT), b_plot(_WIDTH, _HEIGHT) 'for red,green and blue channels
18.  
19. DIM canvas&, curr
20. canvas& = _NEWIMAGE(_WIDTH, _HEIGHT, 32)
21.  
22. RANDOMIZE TIMER
23.  
24. PRINT "Enter the maximum iteration for each pixel : "
25. INPUT temp
26. DIM SHARED seq(temp) AS vec2
27.  
28. CLS
29.  
30. rf! = RND: gf! = RND: bf! = RND
31.  
32.  
33. max_itr = 1000000
34. cx1 = -2: cy1 = -2
35. cx2 = 2: cy2 = 2
36. inc = 1
37.  
38. 'for speed
39. '$CHECKING:OFF
40. DO
41.  
42.     WHILE itr < max_itr
43.         FOR b = 0 TO 15
44.             x = RND * _WIDTH
45.             y = RND * _HEIGHT
46.  
47.             xx = map(x, 0, _WIDTH, cx1, cx2)
48.             yy = map(y, 0, _HEIGHT, cy1, cy2)
49.  
50.             n = num_itr(xx, yy, temp) 'only those points, which escape to infinity are useful.
51.  
52.             'we are now going to store density of each series of this point.
53.             'I've used 3 arrays for color and in each array, the density is stored for different number of iterations.
54.             IF n > 1 THEN
55.                 FOR j = 0 TO n - 2
56.                     px = ABS(INT(map(seq(j).x, cx1, cx2, 1, _WIDTH - 1)))
57.                     py = ABS(INT(map(seq(j).y, cy1, cy2, 1, _HEIGHT - 1)))
58.                     IF px > 0 AND px < _WIDTH AND py > 0 AND py < _HEIGHT THEN
59.  
60.                         IF n < temp * rf! THEN r_plot(px, py) = r_plot(px, py) + inc
61.                         IF n < temp * bf! THEN b_plot(px, py) = b_plot(px, py) + inc
62.                         IF n < temp * gf! THEN g_plot(px, py) = g_plot(px, py) + inc
63.                         PSET (py, px), _RGB(0, 0, 0)
64.                         PSET (py, px), _RGB(r_plot(px, py), g_plot(px, py), b_plot(px, py))
65.                     END IF
66.                 NEXT
67.             END IF
68.         NEXT
69.         itr = itr + 1
70.         IF itr MOD 10000 = 0 THEN _DISPLAY
71.         rendered = 1
72.     WEND
73.     IF rendered = 1 THEN
74.         _FREEIMAGE canvas&
75.         canvas& = _COPYIMAGE(0)
76.         rendered = 0
77.     END IF
78.  
79.     WHILE _MOUSEINPUT: WEND
80.     CLS
81.     _PUTIMAGE , canvas&
82.  
83.     bx1 = _MOUSEX - 50: by1 = _MOUSEY - 50
84.     bx2 = _MOUSEX + 50: by2 = _MOUSEY + 50
85.  
86.     LINE (bx1, by1)-(bx2, by2), _RGB(255, 255, 0), B
87.  
88.     IF _MOUSEBUTTON(1) THEN
89.  
90.         itr = 0
91.         bx1 = -2 * (bx1 / _WIDTH): bx2 = ABS(bx1)
92.         by1 = bx1: by2 = bx2
93.  
94.         cx1 = bx1: cx2 = bx2: cy1 = by1: cy2 = by2
95.         ERASE r_plot, g_plot, b_plot
96.         DIM r_plot(_WIDTH, _HEIGHT), g_plot(_WIDTH, _HEIGHT), b_plot(_WIDTH, _HEIGHT)
97.         CLS
98.         inc = inc + .1
99.     END IF
100.  
101.  
102.     _DISPLAY
103.     _LIMIT 40
104. LOOP
105.  
106. FUNCTION num_itr (x_0, y_0, max_itr)
107.     x = 0: y = 0
108.     x_n_1 = 0: y_n_1 = 0
109.  
110.     FOR i = 0 TO max_itr
111.         x_n_1 = x * x - y * y + x_0
112.         y_n_1 = 2 * x * y + y_0
113.         x = x_n_1
114.         y = y_n_1
115.         seq(i).x = x: seq(i).y = y
116.         IF ABS(x * x + y * y) > 4 THEN num_itr = i: EXIT FUNCTION
117.     NEXT
118. END FUNCTION
119.  
120. FUNCTION map! (value!, minRange!, maxRange!, newMinRange!, newMaxRange!)
121.     map! = ((value! - minRange!) / (maxRange! - minRange!)) * (newMaxRange! - newMinRange!) + newMinRange!
122. END FUNCTION
123.  
124.  

Selecting the area for zooming.


 

Zoomed Fractal

[_vince]

Nice job, hashish. With this fractal, unlike mandelbrot and other escape time fractals, the pixels in your current field of zoom are not the only ones contributing to it.  So you have to guess with some math tricks which pixels are contributing to the current field of vision otherwise you lose density and details.  There's a program that seems to do it well and is GPU accelerated:

https://fractalforums.org/other/55/buddhabrot-magnifier-a-realtime-buddhabrot-zoomer/384

[Ashish]

@Vincent
Hi. I tried to store the the pixels which are contributing to other pixel in a 3-dimensional array (which was actually cumbersome) for adding
more detail as we zoom. However, it broke, instead of working. Maybe, I will try to make a good one with Advance OpenGL and OpenCL in future
just like the GPU program you mention.