Interactive Mandelbrot Set using the Classical Iteration Method

Luis Villaseñor

lvillasen@gmail.com

Usage

Use the left buttom to draw a square to zoom into

Point and click with the right buttom to magnify by a factor of 10

Click with the left button on the rigth side of the image to randomly change the colormap

Click with right button on the right side of the image to set the default colormap

Click on the left side of the image to restart with the full Mandelbrot set

Press the up/down arrow to increase/decrease the maximum number of iterations

Press the right/left arrow to increase/decrease the number of pixels

Type a number from 1-9 to set power index of the iteration formula

Type 'f' to toggle full-screen mode

Type 's' to save the image

License of this example:

GPLv3

Date:

2016-02-09

PyCUDA version:

2015.1.3

   1 # Interactive Mandelbrot Set Accelerated using PyCUDA
   2 # Classical Iteration Method
   3 # Luis Villasenor
   4 # lvillasen@gmail.com
   5 # 2/8/2016
   6 # Licence: GPLv3
   7 
   8 # Usage
   9 # Use the left buttom to draw a square to zoom into 
  10 
  11 # Point and click with the right buttom to magnify by a factor of 10
  12 
  13 # Click with the left button on the rigth side of the 
  14 # image to randomly change the colormap
  15 
  16 # Click with right button on the right side of the image to set the default colormap
  17 
  18 # Click on the left side of the image to restart with the full Mandelbrot set
  19 
  20 # Press the up/down arrow to increase/decrease the maximum number of iterations
  21 
  22 # Press the right/left arrow to increase/decrease the number of pixels
  23 
  24 # Type a number from 1-9 to set power index of the iteration formula
  25 
  26 # Type 'f' to toggle full-screen mode
  27 
  28 # Type 's' to save the image
  29 
  30 
  31 import pycuda.driver as drv
  32 import pycuda.tools
  33 import pycuda.autoinit
  34 import numpy as np
  35 from pylab import cm as cm
  36 import matplotlib.pyplot as plt
  37 from matplotlib.widgets import RectangleSelector
  38 from matplotlib.patches import Rectangle
  39 from pycuda.compiler import SourceModule
  40 global N,n_block,n_grid,x0,y0,side,L,M,power
  41 L=400;
  42 N=800;n_block=16;n_grid=int(N/16);
  43 N=n_block*n_grid;
  44 x0=-.5;y0=0.
  45 side=3.0
  46 i_cmap=49
  47 power=2
  48 fig = plt.figure(figsize=(12,12))
  49 fig.suptitle('Interactive Mandelbrot Set, Accelerated with PyCUDA')
  50 ax = fig.add_subplot(111)
  51 cmaps=[m for m in cm.datad if not m.endswith("_r")]
  52 N,x0,y0,side,L,power
  53 
  54 mod = SourceModule("""
  55     #include <stdio.h>
  56     #include <pycuda-complex.hpp>
  57     #include <math.h>
  58     typedef   pycuda::complex<double> pyComplex;
  59 __device__ float norma(pyComplex z){
  60     return norm(z);
  61 }
  62 __global__ void mandelbrot(double x0, double y0,double side, int L,int power,int *M)
  63 {
  64     int n_x = blockDim.x*gridDim.x;
  65     int idx = threadIdx.x + blockDim.x*blockIdx.x;
  66     int idy = threadIdx.y + blockDim.y*blockIdx.y;
  67     int threadId = idy*n_x+idx;
  68     float delta = side/n_x;
  69     pyComplex c( x0-side/2.+delta*idx,y0-side/2.+delta*idy);
  70     pyComplex z( x0-side/2.+delta*idx,y0-side/2.+delta*idy);
  71     int h = 0;
  72     float R = 2.0;
  73     while( h<L && norma(z)<R){
  74         z=pow(z,power)+c;
  75         h+=1;
  76     }
  77     M[threadId]=h;
  78 }
  79 """)
  80 M = np.zeros((N,N)).astype(np.int32)
  81 func = mod.get_function("mandelbrot")
  82 func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))
  83 def zoom_on_square(eclick, erelease):
  84     'eclick and erelease are the press and release events'
  85     global N,side,x0,y0,myobj,M,power
  86     x1, y1 = min(eclick.xdata,erelease.xdata),min( eclick.ydata,erelease.ydata)
  87     x2, y2 = max(eclick.xdata,erelease.xdata),max( eclick.ydata,erelease.ydata)
  88     #print(" The button you used were: %s %s" % (eclick.button, erelease.button))
  89     #print ' Nx=%d, Ny=%d, x0=%f, y0=%f'%(x1, y1, x0,y0)
  90     #print ' Nx=%d, Ny=%d, x0=%f, y0=%f'%(x2, y2, x0,y0)
  91     x_1=x0+side*(x1-N/2.)/N
  92     y_1=y0+side*(y1-N/2.)/N
  93     x_2=x0+side*(x2-N/2.)/N
  94     y_2=y0+side*(y2-N/2.)/N
  95     x0=(x_2+x_1)/2.
  96     y0=(y_2+y_1)/2.
  97     side=side*(x2-x1+y2-y1)/N/2 # Average of the 2 rectangle sides
  98     func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))    
  99     myobj = plt.imshow(M,origin='lower',cmap=cmaps[i_cmap])
 100     myobj.set_data(M)
 101     ax.add_patch(Rectangle((1 - .1, 1 - .1), 0.2, 0.2,alpha=1, facecolor='none',fill=None, ))
 102     ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 103     plt.draw()
 104 
 105 def key_selector(event):
 106     global N,side,x0,y0,myobj,M,power,L,i_cmap,n_grid
 107     #print(' Key pressed.')
 108     if event.key == u'up':  # Increase max number of iterations
 109         L=int(L*1.2);
 110         print("Maximum number of iterations changed to %d" % L)
 111         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 112         myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 113         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 114         plt.draw()
 115     if event.key == u'down':  # Decrease max number of iterations
 116         L=int(L/1.2);
 117         print("Maximum number of iterations changed to %d" % L)
 118         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 119         myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 120         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 121         plt.draw()
 122     if event.key == u'right':  # Increase  number of pixels
 123         N=int(N*1.2);
 124         n_grid=int(N/16.);
 125         N=n_block*n_grid;
 126         M = np.zeros((N,N)).astype(np.int32)
 127         print("Number of pixels per dimension changed to %d" % N)
 128         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 129         myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 130         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 131         plt.draw()
 132     if event.key == u'left':  # Decrease  number of pixels
 133         N=int(N/1.2);
 134         n_grid=int(N/16.);
 135         N=n_block*n_grid;
 136         M = np.zeros((N,N)).astype(np.int32)
 137         print("Number of pixels per dimension changed to %d" % N)
 138         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 139         myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 140         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 141         plt.draw()
 142     if event.key in ['1','2','3','4','5','6','7','8','9'] :  # Decrease  number of pixels
 143         power=int(event.key)
 144         if power <10 and power >0 : 
 145             print("Power index set to %d" % power)
 146             i_cmap=49
 147             side=3.0; x0=-.5;y0=0.;L=200;
 148             func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))            
 149             myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 150             ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 151             plt.draw()
 152 
 153 key_selector.RS = RectangleSelector(ax, zoom_on_square,
 154                                        drawtype='box', useblit=True,
 155                                        button=[1, 3],  # don't use middle button
 156                                        minspanx=5, minspany=5,
 157                                        spancoords='pixels')
 158 #                                     interactive=False)
 159 def zoom_on_point(event):
 160     global N,side,x0,y0,myobj,L,M,i_cmap,power
 161     #print(" Button pressed: %d" % (event.button))
 162     #print(' event.x= %f, event.y= %f '%(event.x,event.y))
 163     if event.button==3 and event.inaxes: # Zoom on clicked point; new side=10% of old side
 164         x1, y1 = event.xdata, event.ydata
 165         x0=x0+side*(x1-N/2.)/N
 166         y0=y0+side*(y1-N/2.)/N
 167         side=side*.1
 168         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 169         myobj = plt.imshow(M,origin='lower',cmap=cmaps[i_cmap])
 170         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 171         plt.draw()
 172     if not event.inaxes and event.x<.3*N : # Click on left side of image to reset to full fractal
 173         power=2; side=3.0; x0=-.5;y0=0.;i_cmap=49
 174         func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))        
 175         myobj = plt.imshow(M,cmap=cmaps[i_cmap],origin='lower')
 176         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 177         plt.draw()
 178     if event.button==1 and not event.inaxes and event.x>.7*N : # Left click on right side of image to set a random colormap
 179         i_cmap_current=i_cmap
 180         i_cmap=np.random.randint(len(cmaps))
 181         if i_cmap==i_cmap_current:
 182             i_cmap-=1
 183             if i_cmap< 0 : i_cmap=len(cmaps)-1
 184         #print("color=",i_cmap) 
 185         myobj = plt.imshow(M,origin='lower',cmap=cmaps[i_cmap])
 186         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 187         plt.draw() 
 188     if event.button==3 and not event.inaxes and event.x>.7*N : # Right click on right side to set default mapolormap 
 189         i_cmap=49
 190         myobj = plt.imshow(M,origin='lower',cmap=cmaps[i_cmap])
 191         ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 192         plt.draw()   
 193 fig.canvas.mpl_connect('button_press_event', zoom_on_point)
 194 fig.canvas.mpl_connect('key_press_event', key_selector)
 195 func(np.float64(x0),np.float64(y0),np.float64(side), np.int32(L),np.int32(power),drv.Out(M),block=(n_block,n_block,1),grid=(n_grid,n_grid,1))
 196 ax.set_title('Side=%.2e, x=%.2e, y=%.2e, %s, L=%d'%(side,x0,y0,cmaps[i_cmap],L))
 197 plt.imshow(M,origin='lower',cmap=cmaps[i_cmap])
 198 plt.show()


CategoryPyCuda CategoryPyCuda

PyCuda/Examples/MandelbrotInteractive (last edited 2016-02-12 16:44:40 by dsl-187-132-143-20-dyn)