Most magic is the Click & Zoom function, gives you the possibility to keep zooming into interesting area's and generate new images.
Mouse: click to zoom in or out (left/right click), keyboard strokes: Color shuffle: x = Gray, c= Color1, v = color2. Manual color control : u-j, i-k, o-l, [-]
s = Save image, q = QUIT
1 / 3
The accelerated mandel() function has #-marked-out algorithms, check them out, or define your own in the code - HAVE FUN
@cuda.jit(device=True) def mandel(x, y, max): c = complex(y, x) z = 0.0j for i in range(max): #z = z*z*z*c - z*z*c +z*z + c #z = z*z*c + z*c + c #z = (z-c)*(z-1) + z*c #z = z*z*z + c #z = (z-1)*(z-1) + z - c #z = z*z/(z-1) + z/(z-1) + c z = z*z + c #orginal Mandelbrot if (z.real*z.real + z.imag*z.imag) >= 4: return i return 1
NEW : Video fractal !!
Stream video example
NEW: Video Fractal Growth
3D Fractal
Added the option to make a 3D MathPlot of the Fractal. Looks nice, but is not so stunning.
Second option is to export a pointcload txt file (xyzrgb format). This can be imported into a 3D CAD program :
Point Cloud in Rhino3D
Still work to do to create an STL export with color...
Video stream with auto zoom in function - Notebook
{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Fractal Movie Generator #\n","### -file for Fractal calculations- ###\n","\n","JV - 2023\n","\n","'Generates' Fractal movie by using openCV \n","Avoiding to use Matplotlib functions due to slow speed and high memory demand \n","CMAP_BUF : color space : hint: keep iteration same size as color space lenght, ie 512 \n","FRAC: fractal array of floats with values between 0 and 1 \n","FOI : Field of Interest plot, marks points of interest where to zoom in (max unique points)\n","IMAGE_BUF: Image buffer, definition kept out of functions to speed up . \n","\n","Using @Jit NUMBA compiler options for CPU acceleration \n","Hint: check your OpenCV video possibilities, this depends on your installed system (Raspi / Windowsa / Linux etc)\n","All coordinates are used in Y-X (like Row-Column) to avoid messing up your intepretations between graph and arrays :)"]},{"cell_type":"code","execution_count":11,"metadata":{},"outputs":[],"source":["import os \n","import cv2\n","from datetime import datetime\n","import numba \n","from numba import jit\n","import numpy as np\n","import scipy as sc\n","import math\n","from numpy.random import default_rng\n","rng = default_rng()\n","import cmath\n","from numba import cuda\n","cuda.select_device(0)\n","global IMAGENR ; IMAGENR=100\n","global THREADS ; THREADS=720;\n","global Gf , Rf, Bf, Zf, If\n","Gf=2.0;Rf=2.0;Bf=2.0;Zf=0.0;If=0\n","global COLOR_DEPTH\n","COLOR_DEPTH=255"]},{"cell_type":"code","execution_count":12,"metadata":{"tags":[]},"outputs":[],"source":["# Fractal kernal routines for Mandelbrot\n","def calculate_mandelbrot_imageC(Ny,Nx,itmax,center,scl,IMG) : # Callable function\n"," x1=center[1]-Nx/(2*scl); x2=center[1]+Nx/(2*scl)\n"," y1=center[0]-Ny/(2*scl); y2=center[0]+Ny/(2*scl)\n"," \n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp)\n"," grx=math.ceil(Nx/blx); gry=math.ceil(Ny/bly)\n"," blockdim = (blx, bly) ; griddim = (grx,gry)\n"," d_image = cuda.to_device(IMG)\n"," mandel_kernel[griddim, blockdim](x1, x2, y1, y2, d_image, itmax)\n"," IMG = d_image.copy_to_host()\n"," return IMG\n","\n","@cuda.jit(device=True)\n","def mandel(x, y, max):\n"," c = complex(y, x)\n"," z = 0.0j\n"," for i in range(max):\n"," #z = z*z*z*c - z*z*c +z*z + c\n"," #z = z*z*c + z*c + c\n"," #z = (z-c)*(z-1) + z*c\n"," #z = z*z*z + c\n"," #z = (z-1)*(z-1) + z - c\n"," #z = z*z/(z-1) + z/(z-1) + c \n"," z = (z*z) + c #orginal manderbrot \n"," if (z.real*z.real + z.imag*z.imag) >= 4:\n"," return i\n"," return max\n"," \n","@cuda.jit\n","def mandel_kernel(min_x, max_x, min_y, max_y, image, iters):\n"," height = image.shape[0]\n"," width = image.shape[1]\n","\n"," pixel_size_x = (max_x - min_x) / width\n"," pixel_size_y = (max_y - min_y) / height\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," real = min_x + (x * pixel_size_x);\n"," for y in range(startY, height, gridY):\n"," imag = min_y + (y * pixel_size_y);\n"," image[y, x] = mandel(real, imag, iters)/iters;\n"," \n","\n","# Field of interest calculator : edge detector in 5x5 pixel grid\n","def calculate_foiC(M,foi): # Callable function\n"," hy = M.shape[0]\n"," wx = M.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp) # optimum is max 512 threads per Block : gtx1080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly) # calculate grid\n"," blockdim = (blx, bly) ; griddim = (grx,gry) # define block and grid\n"," m_image = cuda.to_device(M)\n"," f_image = cuda.to_device(foi)\n"," foi_kernel[griddim, blockdim](m_image,f_image)\n"," foi=f_image.copy_to_host()\n"," #M=m_image.copy_to_host()\n"," return foi \n","\n","@cuda.jit\n","def foi_kernel(fimagein, fimageout):\n"," height = fimagein.shape[0]\n"," width = fimagein.shape[1]\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," for y in range(startY, height, gridY):\n"," p = fimagein[y,x] ;dev=0\n"," for t in range (-2,3,1) :\n"," for u in range (-2,3,1) :\n"," if y+t>=0 and y+t<height and x+u>=0 and x+u<width :\n"," dev = dev + (p-fimagein[y+t,x+u])*(p-fimagein[y+t,x+u])\n"," fimageout[y, x] = dev/25 # this should be a square root, but for foi indication the sum^2 is enough\n"," \n","\n","# Filter Fractal array\n","def filter_fracC(fr,cv,fmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= fmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080 , 1024 RTX3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," f_map = cuda.to_device(fmap)\n"," filter_kernel[griddim, blockdim](f_image,c_image,f_map)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n"," \n"," \n","@cuda.jit\n","def filter_kernel(frac,cvimg,fmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," fy = fmap.shape[0]\n"," fx = fmap.shape[1]\n"," tt=0;uu=0;\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," p = frac[y,x]; \n"," av=0; wg=0\n"," if p<0.2 or p>0.80:\n"," for t in range (0,fy,1) :\n"," tt = int(t-(fy-1)/2);\n"," for u in range (0,fx,1) :\n"," uu= int(u-(fx-1)/2);\n"," if (y+tt)>=0 and (y+tt)<hy and (x+uu)>=0 and (x+uu)<wx :\n"," wg = wg + fmap[t,u];\n"," av = av + frac[y+tt,x+uu]*fmap[t,u];\n"," cvimg[y,x] = av/wg # this should be a square root, but for foi indication the sum^2 is enough\n"," else:\n"," cvimg[y,x] = p;\n"," \n","\n","\n","# Map fractal array into CV2-image map with R-G-B coding using Color Map\n","def cmap_cvimageC(fr,cv,cmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= cmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," m_image = cuda.to_device(cmap)\n"," cmap_kernel[griddim, blockdim](f_image,c_image,m_image)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n"," #return cv2.cvtColor(np.array(cv), cv2.COLOR_RGB2BGR)\n"," \n","@cuda.jit\n","def cmap_kernel(frac,cvimg,cmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," rr = cmap.shape[0]\n"," cd = cmap.shape[1]\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," index = int(frac[y,x]*(rr-1))\n"," index=index%rr #protection to avoid faul indexes\n"," for z in range(startZ, cd, gridZ):\n"," cvimg[y,x,z]=cmap[index,z]\n","\n"," \n","# calculate zoom point by splitting most interesting point of the FOI into an array, then select point, calculate new center.\n","def calculate_zoompoint(foi,centero,resy,resx, scale, slce) :\n"," foi_r = np.where(foi > np.max(foi)*0.3 ) # results in 2 1 dimensionla arrays col:row\n"," s=np.array([0])\n"," if slce == 0 : s = np.random.uniform(0,foi_r[:][0].size-1,2)\n"," else : s[0]= foi_r[:][0].size*(slce%1)\n"," xindex=foi_r[1][int(s[0])]\n"," yindex=foi_r[0][int(s[0])]\n"," centern= ( centero[0] + (yindex-resy/2)/scale, centero[1] + (xindex-resx/2)/scale, yindex,xindex)\n"," return centern # center results in 4 points: real y/x: -> fractal points and mapped y/x -> indexed points of the array\n","\n","\n","# Complex function to generate 16bit ColorMaps for Fractals by sinus variations. \n","# Seed =0 : Grayscale, Seed = odd : white color Seed = even: black color\n","def create_cmap(buf,seed) :\n"," global Rf; global Bf; \n"," global Gf; global Zf ; global If\n"," r=0;g=0;b=0;inv=0;\n"," x=buf.shape[0]\n"," if seed == 0 : # seed 0 is Gray scale\n"," r=4;g=4;b=4;inv=0;\n"," else : \n"," if seed == 99 : # color frequency by keyboard\n"," #do\n"," r=Rf\n"," g=Gf\n"," b=Bf\n"," inv=If\n"," else :\n"," # other seed is random: \n"," r= np.random.uniform(0.85,1.95,1)\n"," g= np.random.uniform(0.85,1.95,1)\n"," b= np.random.uniform(0.85,1.95,1)\n"," Rf=r;Bf=b;Gf=g;If=inv;\n"," for i in range(0,x,1):\n"," if i*1*math.pi/x < 1*math.pi/6 or i*1*math.pi/x > 5*math.pi/6 : \n"," l1= math.cos(i*3*math.pi/x-math.pi/2) # lightness mask = halve 4*sine, cut off at 1\n"," else : l1= 0.75+ math.cos(i*3*math.pi/x-math.pi/2)/4 # lightness mask = halve 4*sine, cut off at 1\n"," if l1>1 : l1=1\n","\n"," buf[x-i-1][0]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*r*math.pi/x-Zf*math.pi/2))/2 )) # R Channel\n"," buf[x-i-1][1]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*g*math.pi/x-Zf*math.pi/2))/2 )) # G Channel\n"," buf[x-i-1][2]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*b*math.pi/x-Zf*math.pi/2))/2 )) # B Channel\n"," r=r*1.001;g=g*1.001;b=b*1.001\n"," \n"," buf.reshape(x,3)\n"," return buf "]},{"cell_type":"code","execution_count":null,"metadata":{},"outputs":[],"source":["\n"]},{"cell_type":"markdown","metadata":{},"source":["## Start Conditions"]},{"cell_type":"code","execution_count":13,"metadata":{"tags":[]},"outputs":[{"name":"stdout","output_type":"stream","text":["0:00:00.523395\n"]},{"data":{"text/plain":["True"]},"execution_count":13,"metadata":{},"output_type":"execute_result"}],"source":["#Start Conditions\n","iterate_max = 512; center = (-0.25,0, 0,0); scale = 200.0; Ny=520; Nx=960; Zf=0.0; # start parameters for iteration max and centerpoint : Y-X coordinates !!\n","BASE = np.array(np.zeros((Ny,Nx), dtype=np.float64)) # define array global, not in function\n","FOI=BASE.copy()\n","IMAGE_BUF = np.zeros((Ny, Nx,3), dtype = np.uint8)\n","CMAP_BUF = np.zeros((iterate_max,3), dtype = np.uint8)\n","CMAP_BUF = create_cmap(CMAP_BUF,2)\n","FMAP_BUF = np.matrix ([\n"," [0.0, 0.0, 0.2, 0.0, 0.0],\n"," [0.0, 0.2, 0.5, 0.2, 0.0],\n"," [0.2, 0.5, 0.0, 0.5, 0.2],\n"," [0.0, 0.2, 0.5, 0.2, 0.0],\n"," [0.0, 0.0, 0.2, 0.0, 0.0]],dtype=np.float64)\n","\n","t= datetime.now()\n","FRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate mandelbrot\n","NFRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate mandelbrot\n","FRAC =filter_fracC(NFRAC,FRAC,FMAP_BUF)\n","t=datetime.now()-t\n","print(t)\n","C=CMAP_BUF.reshape(1,iterate_max,3)\n","C=np.repeat(C,32,axis=0)\n","cv2.imwrite(\"CMAP.png\", C)\n","frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF)\n","cv2.imwrite(\"MandelOrg.png\", frame) \n"," "]},{"cell_type":"markdown","metadata":{},"source":["## Calculate point of interest, 8 iterations"]},{"cell_type":"code","execution_count":14,"metadata":{},"outputs":[{"data":{"text/plain":["True"]},"execution_count":14,"metadata":{},"output_type":"execute_result"}],"source":["for i in range(0,8,1): # iterations for start - fibnd interesting random area \n"," FOI=calculate_foiC(FRAC,FOI)\n"," center = calculate_zoompoint(FOI,center,Ny,Nx,scale,0)\n"," scale = scale*1.5; # zoom in\n"," FRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,FRAC) # calculate mandelbrot\n","frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF)\n","cv2.imwrite(\"Videostart\"+str(i)+\".png\", frame) \n","#cv2.imshow('frame',frame)"]},{"cell_type":"markdown","metadata":{},"source":["## Stream Video : 7x Panning of 180 zoom steps"]},{"cell_type":"code","execution_count":10,"metadata":{"tags":[]},"outputs":[],"source":["frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF) # first FRAC frame\n","#height, width,c = frame.shape\n","ncenter=center\n","\n","# Define the codec and create VideoWriter object.The output is stored in 'outpy.mp4' file.\n","#out= cv2.VideoWriter('Fractal.mp4', cv2.VideoWriter_fourcc(*'DIVX'), 60, (Nx,Ny))\n","\n","# Define the codec and create VideoWriter object.The output is stored in 'outpy.avi' file.\n","out = cv2.VideoWriter('Fractal.avi',cv2.VideoWriter_fourcc('M','J','P','G'), 40, (Nx,Ny))\n"," \n","for j in range(0,8,1) :\n"," FOI=calculate_foiC(FRAC,FOI)\n"," ncenter = calculate_zoompoint(FOI,center,Ny,Nx,scale,0 )\n"," #print(\"Run \",j,\" - Scale \",scale,\" - Center \",center[0:2])\n"," zstep=140\n"," ystep = (ncenter[0]-center[0])/zstep # panning steps\n"," xstep = (ncenter[1]-center[1])/zstep # panning steps\n"," for i in range(0,zstep,1): # iterations \n"," cuda.synchronize()\n"," scale = scale*1.015; # slowly zoom in\n"," Zf = (Zf+0.005);\n"," CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," center = (center[0] + ystep, center[1] + xstep,0,0); # pan to coordinate\n"," NFRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate mandelbrot\n"," FRAC =filter_fracC(NFRAC,FRAC,FMAP_BUF)\n"," frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF)\n"," cv2.imshow('frame2',frame) \n"," out.write(frame)\n"," key = cv2.waitKey(1) & 0xFF\n"," if key == ord(\"q\"): # if the `q` key was pressed, break from the loop\n"," break \n","out.release()\n","cv2.destroyAllWindows() "]}],"metadata":{"anaconda-cloud":{},"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.3"}},"nbformat":4,"nbformat_minor":4}
CUDAfractal_Click.ipynb
Python
Click and Zoom Notebook
{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Fractal Picture Generator #\n","### python -file for Fractal calculations Cuda accelerated ###\n","\n","JV - 2023/7\n","\n","'Generates' Fractal pictures by using openCV and Cuda\n","Avoiding to use Matplotlib functions due to slow speed and high memory demand \n","Image Matrices / Numpy's :\n","CMAP_BUF : color space \n","FRAC / NFRAC: fractal array of floats64 with relative values between 0 and 1 \n","FOI : Field of Interest plot, marks points of interest where to zoom in (max unique points)\n","IMAGE_BUF: Image buffer, definition kept out of functions to speed up . \n","\n","Using @Jit NUMBA compiler options for CPU acceleration - install the cuda tookit 'conda install cudatoolkit' \n","Check your OpenCV video possibilities, this depends on your installed system (Raspi / Windowsa / Linux etc)\n","All coordinates are used in Y-X (like Row-Column) to avoid messing up your intepretations between graph and arrays :)\n","\n","Click and zoom functions:\n","Mouse: click to zoom in or out (left/right click).\n","* Color shuffle: x = Grayscale, c = New Color map , v = inverse color space, f = filer on/of\n","* Manual color panning: red u-j, green i-k, blue o-l \n","* \\[-\\] shift colors\n","* s = save image\n","* r = reset color space\n","* q = QUIT\n","\n","See def mandel(x, y, max) function to choose different complex functions"]},{"cell_type":"markdown","metadata":{},"source":["## Load Libraries and global variables"]},{"cell_type":"code","execution_count":18,"metadata":{"tags":[]},"outputs":[],"source":["import os \n","import cv2\n","from datetime import datetime\n","import numba \n","from numba import jit\n","import numpy as np\n","import scipy as sc\n","import math\n","from numpy.random import default_rng\n","rng = default_rng()\n","import cmath\n","from numba import cuda\n","cuda.select_device(0)\n","global IMAGENR; IMAGENR=0\n","global THREADS; THREADS=800 # RTX3080 max threads\n","global CLICKL\n","CLICKL=np.array(np.zeros((1,2), dtype=np.int32))\n","global CLICKR\n","CLICKR=np.array(np.zeros((1,2), dtype=np.int32))\n","#global variable for color map generator\n","global Gf, Rf, Bf\n","global Zf, If, Ff\n","Gf=2.0;Rf=2.0;Bf=2.0;Zf=2.0;If=0; Ff=0; # Color spacve parameters\n","global COLOR_DEPTH \n","COLOR_DEPTH = 255; # using 8 bits per color"]},{"cell_type":"markdown","metadata":{},"source":["## Function Declarations"]},{"cell_type":"code","execution_count":19,"metadata":{"tags":[]},"outputs":[],"source":["# Fractal kernal routines for Mandelbrot\n","def calculate_mandelbrot_imageC(Ny,Nx,itmax,center,scl,IMG) : \n"," x1=center[1]-Nx/(2*scl); x2=center[1]+Nx/(2*scl)\n"," y1=center[0]-Ny/(2*scl); y2=center[0]+Ny/(2*scl)\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp) # prepare Cuda RTX3080: warp = 32, Threads = 1024\n"," grx=math.ceil(Nx/blx); gry=math.ceil(Ny/bly)\n"," blockdim = (blx, bly) ; griddim = (grx,gry)\n"," d_image = cuda.to_device(IMG) #send FRAC image to Device\n"," mandel_kernel[griddim, blockdim](x1, x2, y1, y2, d_image, itmax) # call cuda function\n"," IMG = d_image.copy_to_host() #copy image back to host\n"," return IMG\n","\n","@cuda.jit(device=True)\n","def mandel(x, y, max):\n"," c = complex(y, x)\n"," z = 0.0j\n"," z1 = 0.0j\n"," z2 = 0.0j\n"," for i in range(max):\n"," #z = z*z1*z2 - z1*z2 +z + c ; z2=z1; z1=z ; # delayed 2 stage \n"," #z = z*z1 - z*c + z1 -c ; z1 = z ; # delayed \n"," #z = z*z1 + c ; z1 = z ; # delayed\n"," #z = z*z*z*z - z*z + c \n"," #z = z*z*z - z*z +z + c\n"," #z = z*z*c + z*c + c\n"," #z = (z-c)*(z-1) + z*c\n"," #z = z*z*z + c # triple comples\n"," #z = (z-1)*(z+1) + c\n"," #z= (z+c)*(z-c) # double mandelbrot\n"," #z = z*z -z +c # buffalo\n"," z = z*z + c # orginal manderbrot \n"," if (z.real*z.real + z.imag*z.imag) >= 4:\n"," return i\n"," return max\n"," \n","@cuda.jit\n","def mandel_kernel(min_x, max_x, min_y, max_y, image, iters):\n"," height = image.shape[0]\n"," width = image.shape[1]\n","\n"," pixel_size_x = (max_x - min_x) / width\n"," pixel_size_y = (max_y - min_y) / height\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," real = min_x + x * pixel_size_x\n"," for y in range(startY, height, gridY):\n"," imag = min_y + y * pixel_size_y \n"," image[y, x] = mandel(real, imag, iters)/iters\n"," \n","\n","# Field of interest calculator : edge detector in 5x5 pixel grid\n","def calculate_foiC(M,foi): # Callable function\n"," hy = M.shape[0]\n"," wx = M.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp) # optimum is max 32 /512 threads per Block : gtx1080, 32 /1024 for RTx 3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly) # calculate grid\n"," blockdim = (blx, bly) ; griddim = (grx,gry) # define block and grid\n"," m_image = cuda.to_device(M)\n"," f_image = cuda.to_device(foi)\n"," foi_kernel[griddim, blockdim](m_image,f_image)\n"," foi=f_image.copy_to_host()\n"," #M=m_image.copy_to_host()\n"," return foi \n","\n","@cuda.jit\n","def foi_kernel(fimagein, fimageout):\n"," height = fimagein.shape[0]\n"," width = fimagein.shape[1]\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," for y in range(startY, height, gridY):\n"," p = fimagein[y,x] ;dev=0\n"," for t in range (-2,3,1) :\n"," for u in range (-2,3,1) :\n"," if y+t>=0 and y+t<height and x+u>=0 and x+u<width :\n"," dev = dev + (p-fimagein[y+t,x+u])*(p-fimagein[y+t,x+u])\n"," fimageout[y, x] = dev/25 # this should be a square root, but for foi indication the sum^2 is enough\n","\n","\n"," \n"," \n","# Map fractal array into CV2-image map with R-G-B coding using Color Map\n","def cmap_cvimageC(fr,cv,cmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= cmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080 , 1024 RTX3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," m_image = cuda.to_device(cmap)\n"," cmap_kernel[griddim, blockdim](f_image,c_image,m_image)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n","\n"," \n","@cuda.jit\n","def cmap_kernel(frac,cvimg,cmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," rr = cmap.shape[0]\n"," cd = cmap.shape[1]\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," index = int(frac[y,x]*(rr-1))\n"," index=index%rr #protection to avoid faul indexes\n"," for z in range(startZ, cd, gridZ):\n"," cvimg[y,x,z]=cmap[index,z]\n","\n","# Filter Fractal array\n","def filter_fracC(fr,cv,fmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= fmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080 , 1024 RTX3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," f_map = cuda.to_device(fmap)\n"," filter_kernel[griddim, blockdim](f_image,c_image,f_map)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n"," \n"," \n","@cuda.jit\n","def filter_kernel(frac,cvimg,fmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," fy = fmap.shape[0]\n"," fx = fmap.shape[1]\n"," tt=0;uu=0;\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," p = frac[y,x]; \n"," av=0; wg=0\n"," if p<0.2 or p>0.80:\n"," for t in range (0,fy,1) :\n"," tt = int(t-(fy-1)/2);\n"," for u in range (0,fx,1) :\n"," uu= int(u-(fx-1)/2);\n"," if (y+tt)>=0 and (y+tt)<hy and (x+uu)>=0 and (x+uu)<wx :\n"," wg = wg + fmap[t,u];\n"," av = av + frac[y+tt,x+uu]*fmap[t,u];\n"," cvimg[y,x] = av/wg # this should be a square root, but for foi indication the sum^2 is enough\n"," else:\n"," cvimg[y,x] = p;\n"," \n"," \n"," \n"," \n","# calculate zoom point by splitting most interesting point of the FOI into an array, then select point, calculate new center.\n","def calculate_zoompoint(foi,centero,resy,resx, scale, slce) :\n"," foi_r = np.where(foi > np.max(foi)/5 ) # results in 2 1 dimensionla arrays col:row\n"," s=np.array([0])\n"," if slce == 0 : s = np.random.uniform(0,foi_r[:][0].size-1,2)\n"," else : s[0]= foi_r[:][0].size*(slce%1)\n"," xindex=foi_r[1][int(s[0])]\n"," yindex=foi_r[0][int(s[0])]\n"," centern= ( centero[0] + (yindex-resy/2)/scale, centero[1] + (xindex-resx/2)/scale, yindex,xindex)\n"," return centern # center results in 4 points: real y/x: -> fractal points and mapped y/x -> indexed points of the array\n","\n","def calculate_zoompointXY(centero,y,x,resy,resx, scale) :\n"," centern= ( centero[0] + (y-resy/2)/scale, centero[1] + (x-resx/2)/scale, y,x)\n"," return centern # center results in 4 points: real y/x: -> fractal points and mapped y/x -> indexed points of the array\n","\n","\n","\n","# Complex function to generate 16bit ColorMaps for Fractals by sinus variations. \n","# Seed =0 : Grayscale, Seed = odd : white color Seed = even: black color\n","def create_cmap(buf,seed) :\n"," global Rf; global Bf; \n"," global Gf; global Zf ; global If\n"," r=0;g=0;b=0;inv=0;\n"," x=buf.shape[0]\n"," if seed == 0 : # seed 0 is Gray scale\n"," r=4;g=4;b=4;inv=0;\n"," else : \n"," if seed == 99 : # color frequency by keyboard\n"," #do\n"," r=Rf\n"," g=Gf\n"," b=Bf\n"," inv=If\n"," else :\n"," # other seed is random: \n"," r= np.random.uniform(0.65,3.5,1)\n"," g= np.random.uniform(0.65,3.5,1)\n"," b= np.random.uniform(0.65,3.5,1)\n"," Rf=r;Bf=b;Gf=g;If=inv;\n"," for i in range(0,x,1):\n"," if i*1*math.pi/x < 1*math.pi/6 or i*1*math.pi/x > 5*math.pi/6 : \n"," l1= math.cos(i*3*math.pi/x-math.pi/2) # lightness mask = halve 4*sine, cut off at 1\n"," else : l1= 0.75+ math.cos(i*3*math.pi/x-math.pi/2)/4 # lightness mask = halve 4*sine, cut off at 1\n"," if l1>1 : l1=1\n","\n"," buf[x-i-1][0]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*r*math.pi/x-Zf*math.pi/2))/2 )) # R Channel\n"," buf[x-i-1][1]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*g*math.pi/x-Zf*math.pi/2))/2 )) # G Channel\n"," buf[x-i-1][2]= abs( COLOR_DEPTH*inv -int (COLOR_DEPTH*l1*(1+math.sin(i*b*math.pi/x-Zf*math.pi/2))/2 )) # B Channel\n"," #buf[x-i-1][3]=1\n"," r=r*1.001;g=g*1.001;b=b*1.001\n"," \n"," buf.reshape(x,3)\n"," return buf \n","\n","# function to display the coordinates of of the points clicked on the image\n","def click_event(event, x, y, flags, params): \n"," global CLICKL\n"," global CLICKR\n"," # checking for left mouse clicks \n"," if event == cv2.EVENT_LBUTTONDOWN: \n"," CLICKL=(y,x)\n"," if event == cv2.EVENT_RBUTTONDOWN: \n"," CLICKR=(y,x)"]},{"cell_type":"markdown","metadata":{},"source":["## Start Conditions"]},{"cell_type":"code","execution_count":20,"metadata":{"tags":[]},"outputs":[{"name":"stdout","output_type":"stream","text":["916244\n"]},{"data":{"text/plain":["True"]},"execution_count":20,"metadata":{},"output_type":"execute_result"}],"source":["#Start Conditions\n","iterate_max = 256; center = (-0.25,0, 0,0); scale = 500.0; Ny=1280; Nx=1920 # start parameters for iteration max and centerpoint : Y-X coordinates !!, Nx and Ny are image resolution\n","Gf=2.3;Rf=2.1;Bf=1.9;Zf=1.0;If=0; Ff=0;\n","FMAP_BUF = np.matrix ([\n"," [0.0, 0.0, 0.1, 0.0, 0.0],\n"," [0.0, 0.2, 0.4, 0.2, 0.0],\n"," [0.1, 0.4, 0.0, 0.4, 0.1],\n"," [0.0, 0.2, 0.4, 0.2, 0.0],\n"," [0.0, 0.0, 0.1, 0.0, 0.0]],dtype=np.float64)\n","\n","BASE = np.array(np.zeros((Ny,Nx), dtype=np.float64)) # define fractal array global, not in function, float64 type\n","FOI=BASE.copy() # define FOI image, same size\n","IMAGE_BUF = np.zeros((Ny, Nx,3), dtype = np.uint8) # define Color image buffer\n","CMAP_BUF = np.zeros((iterate_max,3), dtype = np.uint8) # define Color Map 8 bit type, x 3 (24 bit RGB) - no alfa channel\n","CMAP_BUF = create_cmap(CMAP_BUF,2) # create a color map\n","C=CMAP_BUF.reshape(1,iterate_max,3)\n","C=np.repeat(C,32,axis=0)\n","cv2.imwrite(\"CMAP.png\", C)\n","\n","t= datetime.now()\n","FRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate mandelbrot\n","NFRAC = FRAC.copy();\n","NFRAC=filter_fracC(FRAC,NFRAC,FMAP_BUF)\n","frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n","t=datetime.now()-t\n","print(t.microseconds)\n","frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF)\n","cv2.imwrite(\"MandelOrg.png\", frame) \n"," "]},{"cell_type":"markdown","metadata":{"tags":[]},"source":["## Generate first zoom point by FOI"]},{"cell_type":"code","execution_count":21,"metadata":{"tags":[]},"outputs":[{"data":{"text/plain":["True"]},"execution_count":21,"metadata":{},"output_type":"execute_result"}],"source":["for i in range(0,8,1) :\n"," FRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate mandelbrot\n"," FOI=calculate_foiC(FRAC,FOI)\n"," center = calculate_zoompoint(FOI,center,Ny,Nx,scale,0 )\n"," scale = scale*2; iterate_max = iterate_max+0.25\n","NFRAC=filter_fracC(FRAC,NFRAC,FMAP_BUF)\n","frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n","cv2.imwrite(\"MandelStart.png\", frame) "]},{"cell_type":"markdown","metadata":{},"source":["## Click and Zoom"]},{"cell_type":"code","execution_count":22,"metadata":{"tags":[]},"outputs":[],"source":["#FMAP_BUF, CMAP_BUF and FRAC already defined = mandel set\n","CLICKL=(0,0)\n","CLICKR=(0,0)\n","nr = 0\n","while 1 :\n"," cuda.synchronize()\n"," FRAC = calculate_mandelbrot_imageC(Ny,Nx,iterate_max,center,scale,BASE) # calculate new mandelbrot Set\n"," if Ff ==1 : \n"," NFRAC=filter_fracC(FRAC,NFRAC,FMAP_BUF) ; # filter\n"," frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n"," else : frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF) \n"," cv2.imshow('image',frame)\n"," \n"," #C=CMAP_BUF.reshape(1,iterate_max,3)\n"," #C=np.repeat(C,32,axis=0)\n"," #cv2.imshow('cmap',C)\n"," \n"," key = cv2.waitKey(1) & 0xFF\n"," if key == ord(\"q\"): # if the `q` key was pressed, break from the loop\n"," break \n"," if key == ord(\"s\"): # if the `s pressed, save picture in PNG format\n"," cv2.imwrite(\"MandelZoom\"+str(IMAGENR)+\".png\", frame) \n"," IMAGENR=IMAGENR+1\n"," if key == ord(\"c\"): # if the `c` create new color map black\n"," CMAP_BUF = create_cmap(CMAP_BUF,1)\n"," if key == ord(\"v\"): # if the `v` toggle inversed map\n"," if If == 0: If=1;\n"," else : If=0;\n"," CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"f\"): # if the `v` toggle filter on of\n"," if Ff == 0: Ff=1;\n"," else : Ff=0; \n"," CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"x\"): # if the `x` set grayscale\n"," Gf=Rf;Bf=Rf; CMAP_BUF = create_cmap(CMAP_BUF,0)\n"," if key == ord(\"u\"): # if the `u` pan color map red to higher freq\n"," Rf=Rf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"j\"): # if the `j` pan color map red to lower freq\n"," Rf=Rf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"i\"): # if the `i` pan colormap green to higher freq\n"," Gf=Gf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"k\"): # if the `k` pan colormap green to lower freq\n"," Gf=Gf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"o\"): # if the `o` pan color map blue to higher freq\n"," Bf=Bf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"l\"): # if the `l` pan color map blue to lower freq\n"," Bf=Bf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"[\"): # if the `[` rotate colormap left\n"," Zf=Zf*1.02; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"]\"): # if the `]' srotate color map right\n"," Zf=Zf*0.98; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"r\"): # reset all color parameters\n"," Gf=2.3;Rf=2.1;Bf=1.9;Zf=1.0;If=0; Ff=0; CMAP_BUF = create_cmap(CMAP_BUF,99) \n"," \n"," cv2.setMouseCallback('image', click_event) \n"," if CLICKL != (0,0) :\n"," center = calculate_zoompointXY(center,CLICKL[0],CLICKL[1],Ny,Nx,scale)\n"," #print(CLICK,center)\n"," scale = scale*3; #iterate_max = iterate_max*1.1 # slowly zoom in\n"," CLICKL=(0,0);\n"," if CLICKR != (0,0) :\n"," center = calculate_zoompointXY(center,CLICKR[0],CLICKR[1],Ny,Nx,scale)\n"," #print(CLICK,center)\n"," scale = scale/3; #iterate_max = iterate_max/1.1 # slowly zoom out\n"," CLICKR=(0,0); \n"," \n","cv2.destroyAllWindows() "]},{"cell_type":"code","execution_count":30,"metadata":{},"outputs":[],"source":["\n"]},{"cell_type":"code","execution_count":48,"metadata":{"tags":[]},"outputs":[],"source":[]},{"cell_type":"code","execution_count":null,"metadata":{"tags":[]},"outputs":[],"source":[]}],"metadata":{"anaconda-cloud":{},"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.3"}},"nbformat":4,"nbformat_minor":4}
CUDAfractal_Grow.ipynb
Python
3rd dimension fractal generator - with video generator
{"cells":[{"cell_type":"markdown","metadata":{},"source":["# Fractal Generator - Growth #\n","### python -file for Fractal calculations Cuda accelerated ###\n","\n","JV - 2023 / 8\n","\n","'Generates' Fractal pictures by using openCV \n","Avoiding to use Matplotlib functions due to slow speed and high memory demand \n","Image Matrices / Numpy's :\n","CMAP_BUF : color space \n","FRAC / NFRAC: fractal array of floats64 with relative values between 0 and 1 \n","FOI : Field of Interest plot, marks points of interest where to zoom in (max unique points)\n","IMAGE_BUF: Image buffer, definition kept out of functions to speed up . \n","\n","Using @Jit NUMBA compiler options for CPU acceleration - install the cuda tookit 'conda install cudatoolkit' \n","Check your OpenCV video possibilities, this depends on your installed system (Raspi / Windowsa / Linux etc)\n","All coordinates are used in Y-X (like Row-Column) to avoid messing up your intepretations between graph and arrays :)\n","\n","Click and zoom functions:\n","Mouse: click to zoom in or out (left/right click).\n","* Color shuffle: x = Grayscale, c = New Color map , v = inverse color space, f = filer on/of\n","* Manual color panning: red u-j, green i-k, blue o-l \n","* \\[-\\] shift colors\n","* a-z PAN 3rd DIMENSION grower < ----------------------------------------\n","* s = save image\n","* r = reset color space\n","* q = QUIT\n","\n","See def mandel(x, y, max) function to choose different complex functions"]},{"cell_type":"markdown","metadata":{},"source":["## Load Libraries and global variables"]},{"cell_type":"code","execution_count":8,"metadata":{"tags":[]},"outputs":[],"source":["import os \n","import cv2\n","from datetime import datetime\n","import numba \n","from numba import jit\n","import numpy as np\n","import scipy as sc\n","import math\n","import time\n","from numpy.random import default_rng\n","rng = default_rng()\n","import cmath\n","from numba import cuda\n","cuda.select_device(0)\n","global IMAGENR; IMAGENR=0\n","global THREADS; THREADS=800 # RTX3080 max threads\n","global CLICKL\n","CLICKL=np.array(np.zeros((1,2), dtype=np.int32))\n","global CLICKR\n","CLICKR=np.array(np.zeros((1,2), dtype=np.int32))\n","#global variable for color map generator\n","global Gf\n","global Rf\n","global Bf\n","global Zf\n","global If\n","global Ff\n","Gf=2.0;Rf=2.0;Bf=2.0;Zf=2.0;If=0; Ff=0; # Color spacve parameters"]},{"cell_type":"markdown","metadata":{},"source":["## Function Declarations"]},{"cell_type":"code","execution_count":9,"metadata":{"tags":[]},"outputs":[],"source":["# Fractal kernal routines for Mandelbrot\n","def calculate_mandelbrot_imageC(Ny,Nx,Zz,itmax,center,scl,IMG) : \n"," x1=center[1]-Nx/(2*scl); x2=center[1]+Nx/(2*scl)\n"," y1=center[0]-Ny/(2*scl); y2=center[0]+Ny/(2*scl) \n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp) # prepare Cuda RTX3080: warp = 32, Threads = 1024\n"," grx=math.ceil(Nx/blx); gry=math.ceil(Ny/bly)\n"," blockdim = (blx, bly) ; griddim = (grx,gry)\n"," d_image = cuda.to_device(IMG) #send FRAC image to Device\n"," mandel_kernel[griddim, blockdim](x1, x2, y1, y2, Zz , d_image, itmax) # call cuda function\n"," IMG = d_image.copy_to_host() #copy image back to host\n"," return IMG\n","\n","@cuda.jit(device=True)\n","def mandel(x, y, u, max):\n"," c = complex(y, x) ; \n"," z= complex(0,0); \n"," \n"," for i in range(max):\n"," z1= z*z + c\n"," z2= z*z - u\n"," \n"," #z = z1*z2 # amoebe 1\n"," z = (z-z1)*(z-z2) + c # amoebe 2 \n"," #z = z*z+c\n"," if ((z.real)**2 + (z.imag)**2) >= 4:\n"," return i\n"," return max\n"," \n","@cuda.jit\n","def mandel_kernel(min_x, max_x, min_y, max_y,max_z, image, iters):\n"," height = image.shape[0]\n"," width = image.shape[1]\n","\n"," pixel_size_x = (max_x - min_x) / width\n"," pixel_size_y = (max_y - min_y) / height\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," real = min_x + x * pixel_size_x\n"," for y in range(startY, height, gridY):\n"," imag = min_y + y * pixel_size_y \n"," image[y, x] = mandel(real, imag, max_z, iters)/iters\n"," \n","\n","# Field of interest calculator : edge detector in 5x5 pixel grid\n","def calculate_foiC(M,foi): # Callable function\n"," hy = M.shape[0]\n"," wx = M.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp) # optimum is max 32 /512 threads per Block : gtx1080, 32 /1024 for RTx 3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly) # calculate grid\n"," blockdim = (blx, bly) ; griddim = (grx,gry) # define block and grid\n"," m_image = cuda.to_device(M)\n"," f_image = cuda.to_device(foi)\n"," foi_kernel[griddim, blockdim](m_image,f_image)\n"," foi=f_image.copy_to_host()\n"," #M=m_image.copy_to_host()\n"," return foi \n","\n","@cuda.jit\n","def foi_kernel(fimagein, fimageout):\n"," height = fimagein.shape[0]\n"," width = fimagein.shape[1]\n","\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n","\n"," for x in range(startX, width, gridX):\n"," for y in range(startY, height, gridY):\n"," p = fimagein[y,x] ;dev=0\n"," for t in range (-2,3,1) :\n"," for u in range (-2,3,1) :\n"," if y+t>=0 and y+t<height and x+u>=0 and x+u<width :\n"," dev = dev + (p-fimagein[y+t,x+u])*(p-fimagein[y+t,x+u])\n"," fimageout[y, x] = dev/25 # this should be a square root, but for foi indication the sum^2 is enough\n","\n","\n"," \n"," \n","# Map fractal array into CV2-image map with R-G-B coding using Color Map\n","def cmap_cvimageC(fr,cv,cmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= cmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080 , 1024 RTX3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," m_image = cuda.to_device(cmap)\n"," cmap_kernel[griddim, blockdim](f_image,c_image,m_image)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n","\n"," \n","@cuda.jit\n","def cmap_kernel(frac,cvimg,cmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," rr = cmap.shape[0]\n"," cd = cmap.shape[1]\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," index = int(frac[y,x]*(rr-1))\n"," index=index%rr #protection to avoid faul indexes\n"," for z in range(startZ, cd, gridZ):\n"," cvimg[y,x,z]=cmap[index,z]\n","\n","# Filter Fractal array\n","def filter_fracC(fr,cv,fmap) : #callable function\n"," hy = fr.shape[0]\n"," wx = fr.shape[1]\n"," cd= fmap.shape[1]\n"," warp = 32\n"," blx = warp; bly = math.ceil(THREADS/warp); blz=1 # optimum is max 512 threads per Block : gtx1080 , 1024 RTX3080\n"," grx=math.ceil(wx/blx); gry=math.ceil(hy/bly); grz= cd # calculate grid\n"," blockdim = (blx, bly, blz) ; griddim = (grx,gry,grz) # define block and grid\n"," f_image = cuda.to_device(fr)\n"," c_image = cuda.to_device(cv)\n"," f_map = cuda.to_device(fmap)\n"," filter_kernel[griddim, blockdim](f_image,c_image,f_map)\n"," #cmap=m_image.copy_to_host()\n"," cv=c_image.copy_to_host()\n"," #fr=f_image.copy_to_host() \n"," return cv\n"," \n"," \n","@cuda.jit\n","def filter_kernel(frac,cvimg,fmap) :\n"," hy = frac.shape[0]\n"," wx = frac.shape[1]\n"," fy = fmap.shape[0]\n"," fx = fmap.shape[1]\n"," tt=0;uu=0;\n"," startX = cuda.blockDim.x * cuda.blockIdx.x + cuda.threadIdx.x\n"," startY = cuda.blockDim.y * cuda.blockIdx.y + cuda.threadIdx.y\n"," startZ = cuda.blockDim.z * cuda.blockIdx.z + cuda.threadIdx.z\n"," gridX = cuda.gridDim.x * cuda.blockDim.x;\n"," gridY = cuda.gridDim.y * cuda.blockDim.y;\n"," gridZ = cuda.gridDim.z * cuda.blockDim.z;\n"," for x in range(startX, wx, gridX):\n"," for y in range(startY, hy, gridY):\n"," p = frac[y,x]; \n"," av=0; wg=0\n"," if p<0.2 or p>0.80:\n"," for t in range (0,fy,1) :\n"," tt = int(t-(fy-1)/2);\n"," for u in range (0,fx,1) :\n"," uu= int(u-(fx-1)/2);\n"," if (y+tt)>=0 and (y+tt)<hy and (x+uu)>=0 and (x+uu)<wx :\n"," wg = wg + fmap[t,u];\n"," av = av + frac[y+tt,x+uu]*fmap[t,u];\n"," cvimg[y,x] = av/wg # this should be a square root, but for foi indication the sum^2 is enough\n"," else:\n"," cvimg[y,x] = p;\n"," \n"," \n"," \n"," \n","# calculate zoom point by splitting most interesting point of the FOI into an array, then select point, calculate new center.\n","def calculate_zoompoint(foi,centero,resy,resx, scale, slce) :\n"," foi_r = np.where(foi > np.max(foi)/5 ) # results in 2 1 dimensionla arrays col:row\n"," s=np.array([0])\n"," if slce == 0 : s = np.random.uniform(0,foi_r[:][0].size-1,2)\n"," else : s[0]= foi_r[:][0].size*(slce%1)\n"," xindex=foi_r[1][int(s[0])]\n"," yindex=foi_r[0][int(s[0])]\n"," centern= ( centero[0] + (yindex-resy/2)/scale, centero[1] + (xindex-resx/2)/scale, yindex,xindex)\n"," return centern # center results in 4 points: real y/x: -> fractal points and mapped y/x -> indexed points of the array\n","\n","def calculate_zoompointXY(centero,y,x,resy,resx, scale) :\n"," centern= ( centero[0] + (y-resy/2)/scale, centero[1] + (x-resx/2)/scale, y,x)\n"," return centern # center results in 4 points: real y/x: -> fractal points and mapped y/x -> indexed points of the array\n","\n","\n","\n","# Complex function to generate 16bit ColorMaps for Fractals by sinus variations. \n","# Seed =0 : Grayscale, Seed = odd : white color Seed = even: black color\n","def create_cmap(buf,seed) :\n"," global Rf; global Bf; \n"," global Gf; global Zf ; global If\n"," r=0;g=0;b=0;inv=0;\n"," color_bits = 256-1; # using 16 bits per color\n"," x=buf.shape[0]\n"," if seed == 0 : # seed 0 is Gray scale\n"," r=4;g=4;b=4;inv=0;\n"," else : \n"," if seed == 99 : # color frequency by keyboard\n"," #do\n"," r=Rf\n"," g=Gf\n"," b=Bf\n"," inv=If\n"," else :\n"," # other seed is random: \n"," r= np.random.uniform(0.65,3.5,1)\n"," g= np.random.uniform(0.65,3.5,1)\n"," b= np.random.uniform(0.65,3.5,1)\n"," Rf=r;Bf=b;Gf=g;If=inv;\n"," for i in range(0,x,1):\n"," if i*1*math.pi/x < 1*math.pi/6 or i*1*math.pi/x > 5*math.pi/6 : \n"," l1= math.cos(i*3*math.pi/x-math.pi/2) # lightness mask = halve 4*sine, cut off at 1\n"," else : l1= 0.75+ math.cos(i*3*math.pi/x-math.pi/2)/4 # lightness mask = halve 4*sine, cut off at 1\n"," if l1>1 : l1=1\n","\n"," buf[x-i-1][0]= abs( color_bits*inv -int (color_bits*l1*(1+math.sin(i*r*math.pi/x-Zf*math.pi/2))/2 )) # R Channel\n"," buf[x-i-1][1]= abs( color_bits*inv -int (color_bits*l1*(1+math.sin(i*g*math.pi/x-Zf*math.pi/2))/2 )) # G Channel\n"," buf[x-i-1][2]= abs( color_bits*inv -int (color_bits*l1*(1+math.sin(i*b*math.pi/x-Zf*math.pi/2))/2 )) # B Channel\n"," #buf[x-i-1][3]=1\n"," r=r*1.001;g=g*1.001;b=b*1.001\n"," \n"," buf.reshape(x,3)\n"," return buf \n","\n","# function to display the coordinates of of the points clicked on the image\n","def click_event(event, x, y, flags, params): \n"," global CLICKL\n"," global CLICKR\n"," # checking for left mouse clicks \n"," if event == cv2.EVENT_LBUTTONDOWN: \n"," CLICKL=(y,x)\n"," if event == cv2.EVENT_RBUTTONDOWN: \n"," CLICKR=(y,x)"]},{"cell_type":"markdown","metadata":{},"source":["## Start Conditions"]},{"cell_type":"code","execution_count":10,"metadata":{"tags":[]},"outputs":[{"name":"stdout","output_type":"stream","text":["736272\n"]},{"data":{"text/plain":["True"]},"execution_count":10,"metadata":{},"output_type":"execute_result"}],"source":["#Start Conditions\n","iterate_max = 256 ; center = (-0.25,0, 0,0); scale = 200.0; Ny=1080; Nx=1080; Nz=480;Zz=-1.06 # start parameters for iteration max and centerpoint : Y-X coordinates !!, Nx and Ny are image resolution\n","Gf=2.3;Rf=2.1;Bf=1.9;Zf=1.0;If=0; Ff=0;\n","FMAP_BUF = np.matrix ([\n"," [0.0, 0.0, 0.1, 0.0, 0.0],\n"," [0.0, 0.2, 0.4, 0.2, 0.0],\n"," [0.1, 0.4, 0.0, 0.4, 0.1],\n"," [0.0, 0.2, 0.4, 0.2, 0.0],\n"," [0.0, 0.0, 0.1, 0.0, 0.0]],dtype=np.float64)\n","\n","MATRIX = np.array(np.zeros((Ny,Nx,Nz), dtype=np.float64)) \n","BASE = np.array(np.zeros((Ny,Nx), dtype=np.float64)) # define fractal array global, not in function, float64 type\n","FOI = np.array(np.zeros((Ny,Nx), dtype=np.float64)) # define FOI image, same size\n","IMAGE_BUF = np.zeros((Ny, Nx,3), dtype = np.uint8) # define Color image buffer\n","CMAP_BUF = np.zeros((iterate_max,3), dtype = np.uint8) # define Color Map 16 bit type, x 3 (RGB) - no alfa channel\n","CMAP_BUF = create_cmap(CMAP_BUF,2) # create a color map\n","C=CMAP_BUF.reshape(1,iterate_max,3)\n","C=np.repeat(C,32,axis=0)\n","cv2.imwrite(\"CMAP.png\", C)\n","\n","t= datetime.now()\n","FRAC = calculate_mandelbrot_imageC(Ny,Nx,Zz,iterate_max,center,scale,BASE) # calculate mandelbrot\n","NFRAC = FRAC.copy();\n","NFRAC=filter_fracC(FRAC,NFRAC,FMAP_BUF)\n","frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n","t=datetime.now()-t\n","print(t.microseconds)\n","frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF)\n","cv2.imwrite(\"MandelOrg.png\", frame) \n"," "]},{"cell_type":"markdown","metadata":{},"source":["## Click and Zoom Fractal Growth"]},{"cell_type":"code","execution_count":11,"metadata":{"tags":[]},"outputs":[],"source":["#FMAP_BUF, CMAP_BUF and FRAC already defined = mandel set\n","CLICKL=(0,0)\n","CLICKR=(0,0)\n","nr = 0\n","while 1 :\n"," cuda.synchronize()\n"," FRAC = calculate_mandelbrot_imageC(Ny,Nx,Zz,iterate_max,center,scale,BASE) # calculate new mandelbrot Set\n"," if Ff ==1 : \n"," NFRAC=filter_fracC(FRAC,NFRAC,FMAP_BUF) ; # filter\n"," frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n"," else : frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF) \n"," cv2.imshow('image',frame)\n"," \n"," #C=CMAP_BUF.reshape(1,iterate_max,3)\n"," #C=np.repeat(C,32,axis=0)\n"," #cv2.imshow('cmap',C)\n"," \n"," key = cv2.waitKey(1) & 0xFF\n"," if key == ord(\"q\"): # if the `q` key was pressed, break from the loop\n"," break \n"," if key == ord(\"s\"): # if the `s pressed, save picture in PNG format\n"," cv2.imwrite(\"MandelZoom\"+str(IMAGENR)+\".png\", frame) \n"," IMAGENR=IMAGENR+1\n"," if key == ord(\"z\"): # if the `z` Pan 3rd dimension < ----------------------------------------------------------------- !!!\n"," Zz = Zz + 1/scale \n"," if key == ord(\"a\"): # if the `a` Pan 3rd dimension < --------------------------------------------------------------- !!!\n"," Zz = Zz - 1/scale \n"," if key == ord(\"c\"): # if the `c` create new color map black\n"," CMAP_BUF = create_cmap(CMAP_BUF,1)\n"," if key == ord(\"v\"): # if the `v` toggle inversed map\n"," if If == 0: If=1;\n"," else : If=0;\n"," CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"f\"): # if the `v` toggle filter on of\n"," if Ff == 0: Ff=1;\n"," else : Ff=0; \n"," CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"x\"): # if the `x` set grayscale\n"," Gf=Rf;Bf=Rf; CMAP_BUF = create_cmap(CMAP_BUF,0)\n"," if key == ord(\"u\"): # if the `u` pan color map red to higher freq\n"," Rf=Rf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"j\"): # if the `j` pan color map red to lower freq\n"," Rf=Rf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"i\"): # if the `i` pan colormap green to higher freq\n"," Gf=Gf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"k\"): # if the `k` pan colormap green to lower freq\n"," Gf=Gf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"o\"): # if the `o` pan color map blue to higher freq\n"," Bf=Bf*1.05; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"l\"): # if the `l` pan color map blue to lower freq\n"," Bf=Bf*0.95; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"[\"): # if the `[` rotate colormap left\n"," Zf=Zf*1.02; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"]\"): # if the `]' srotate color map right\n"," Zf=Zf*0.98; CMAP_BUF = create_cmap(CMAP_BUF,99)\n"," if key == ord(\"r\"): # reset all color parameters\n"," Gf=2.3;Rf=2.1;Bf=1.9;Zf=1.0;If=0; Ff=0; CMAP_BUF = create_cmap(CMAP_BUF,99) \n"," \n"," cv2.setMouseCallback('image', click_event) \n"," if CLICKL != (0,0) :\n"," center = calculate_zoompointXY(center,CLICKL[0],CLICKL[1],Ny,Nx,scale)\n"," #print(CLICK,center)\n"," scale = scale*3; #iterate_max = iterate_max*1.1 # slowly zoom in\n"," CLICKL=(0,0);\n"," if CLICKR != (0,0) :\n"," center = calculate_zoompointXY(center,CLICKR[0],CLICKR[1],Ny,Nx,scale)\n"," #print(CLICK,center)\n"," scale = scale/3; #iterate_max = iterate_max/1.1 # slowly zoom out\n"," CLICKR=(0,0); \n"," \n","cv2.destroyAllWindows() "]},{"cell_type":"markdown","metadata":{},"source":["## EXTRA : Fractal Growth Video"]},{"cell_type":"code","execution_count":12,"metadata":{"tags":[]},"outputs":[],"source":["frame=cmap_cvimageC(FRAC,IMAGE_BUF,CMAP_BUF) # first FRAC frame\n","\n","# Define the codec and create VideoWriter object.The output is stored in 'outpy.mp4' file.\n","#out= cv2.VideoWriter('Fractal3D.mp4', cv2.VideoWriter_fourcc(*'DIVX'), 60, (Nx,Ny))\n","out = cv2.VideoWriter('Fractal3D.avi',cv2.VideoWriter_fourcc('M','J','P','G'), 40, (Nx,Ny))\n","#out = cv2.VideoWriter('fraxvid.avi', cv2.VideoWriter_fourcc('X','V','I','D'), 40, (Nx,Ny))\n","\n","Zo=Zz \n","for j in range(0,Nz,1) :\n"," cuda.synchronize()\n"," Zz= Zo + (j-Nz/2)/(scale)\n"," NFRAC = calculate_mandelbrot_imageC(Ny,Nx,Zz,iterate_max,center,scale,BASE) # calculate mandelbrot\n"," #FRAC =filter_fracC(NFRAC,FRAC,FMAP_BUF)\n"," MATRIX[:,:,j]=NFRAC # map frame into Matrix\n"," frame=cmap_cvimageC(NFRAC,IMAGE_BUF,CMAP_BUF)\n"," cv2.imshow('frame2',frame) \n"," out.write(frame)\n"," key = cv2.waitKey(1) & 0xFF\n"," if key == ord(\"q\"): # if the `q` key was pressed, break from the loop\n"," break \n","out.release()\n","cv2.destroyAllWindows() "]},{"cell_type":"markdown","metadata":{},"source":["## WRITE Frame to 3D pointcloud XYZrgb TXT-FILE ##"]},{"cell_type":"code","execution_count":27,"metadata":{},"outputs":[],"source":["l=[]\n","colorbits = 256-1 # color bits of the CMAP matrix : 16 bits\n","rr = CMAP_BUF.shape[0]\n","for i in range(0,Ny,1):\n"," for j in range(0,Nx,1) :\n"," fr= FRAC[i, j]\n"," index = fr*(rr-1) ;\n"," index = int(index%rr) # color index \n"," red = int(CMAP_BUF[index,0]*255/colorbits)\n"," green = int(CMAP_BUF[index,1]*255/colorbits)\n"," blue = int(CMAP_BUF[index,2]*255/colorbits)\n"," if i==Ny-1 and j==Nx-1 : l.append('%d,%d,%d,%d,%d,%d' %(i, j, 100*fr, red , green, blue ) )\n"," else : l.append('%d,%d,%d,%d,%d,%d\\n' %(i, j, 100*fr, red , green, blue ) )\n","\n","with open('fractal.xyz', 'w') as f:\n"," f.write(\"\".join(l))\n","f.close()\n","\n"]}],"metadata":{"anaconda-cloud":{},"kernelspec":{"display_name":"Python 3 (ipykernel)","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.11.3"}},"nbformat":4,"nbformat_minor":4}
_KhLqgQRnE8.gif?auto=format%2Ccompress&gifq=35&w=400&h=300&fit=min)
Comments