SHARE
TWEET

Mandelbrot set explorer

a guest Feb 16th, 2013 91 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. #include <iostream>
  2. #include "timer.h"
  3. #include <stdio.h>
  4. #include <opencv2/core/core.hpp>
  5. #include <opencv2/highgui/highgui.hpp>
  6. #include <opencv2/opencv.hpp>
  7. #include <cuda.h>
  8. #include <cuda_runtime.h>
  9. #include <string>
  10.  
  11. using namespace std;
  12.  
  13. /*
  14. * Global variables
  15. */
  16. int numCols = 800, numRows = 600; // Image size
  17. double cx = 0, cy = 0, zoom = 100.0f; // Center of the image and zoom
  18. uchar *h_image; // The image
  19. int iterationMax = 100; // Maximum number of iterations in the Mandelbrot loop
  20. int framesPerZoom = 10; // Number of frames in the zoom animation
  21.  
  22.  
  23.  
  24. /*
  25. * Operations on complex numbers
  26. */
  27. __device__ void cMult(double *ar, double *ai, double *br, double *bi, double  *cr, double  *ci)
  28. {
  29.         *cr = *ar*(*br) - *ai*(*bi);
  30.         *ci = *ar*(*bi) + *ai*(*br);
  31. }
  32. __device__ void cAdd(double *ar, double *ai, double *br, double *bi, double  *cr, double  *ci)
  33. {
  34.         *cr = *ar + *br;
  35.         *ci = *ai + *bi;
  36. }
  37. __device__ void cMin(double *ar, double *ai, double *br, double *bi, double  *cr, double  *ci)
  38. {
  39.         *cr = *ar - *br;
  40.         *ci = *ai - *bi;
  41. }
  42.  
  43. /*
  44. * Fractal functions
  45. */
  46. // Z^2 + C
  47. __device__ void fractal1(double *x, double *y, double a, double  b)
  48. {
  49.         double auxX, auxY;
  50.         cMult(x, y, x, y, &auxX, &auxY);
  51.         *x = auxX + a;
  52.         *y = auxY + b;
  53. }
  54. // Z^4 + C
  55. __device__ void fractal2(double *x, double *y, double  a, double  b)
  56. {
  57.         double auxX, auxY;
  58.         cMult(x, y, x, y, &auxX, &auxY);
  59.         *x = auxX;
  60.         *y = auxY;
  61.         cMult(x, y, x, y, &auxX, &auxY);
  62.         *x = auxX + a;
  63.         *y = auxY + b;
  64. }
  65.  
  66.  
  67. /*
  68. * Kernel function
  69. */
  70. __global__ void mandelbrot(uchar* const d_image, int numRows, int numCols, double cx, double cy, double zoom, int iterationMax)
  71. {
  72.         // Init
  73.         const int2 thread_2D_pos = make_int2( blockIdx.x * blockDim.x + threadIdx.x,
  74.                 blockIdx.y * blockDim.y + threadIdx.y);
  75.  
  76.         const int thread_1D_pos = (thread_2D_pos.y * numCols + thread_2D_pos.x)*3;
  77.  
  78.         // Exit if we are out of the image
  79.         if(thread_2D_pos.x >= numCols || thread_2D_pos.y >= numRows) return;
  80.  
  81.         // Mandelbrot loop
  82.         /*
  83.                 U0 = 0
  84.                 Un+1 = Un*Un + c
  85.                 where c = a + b*i
  86.                 We stop when |Un|>2 or after iterationMax iterations
  87.  
  88.                 You can find many similar formulas here : http://www.lifesmith.com/formulas.html
  89.         */
  90.         int i = 0;
  91.         double a = cx + (thread_2D_pos.x-numCols/2)/zoom, b = cy + (thread_2D_pos.y- numRows/2)/zoom, x = 0, y = 0;
  92.         for(i = 0 ; i < iterationMax ; i++)
  93.         {
  94.                 fractal1(&x, &y, a, b);
  95.  
  96.                 if(x*x + y*y > 4) break;
  97.         }
  98.  
  99.         // Color conversion
  100.         i = (1.*i)/iterationMax*255.0;
  101.         if(i/64 == 0)
  102.         {
  103.                 d_image[thread_1D_pos] = 0;
  104.                 d_image[thread_1D_pos+1] = i*4;
  105.                 d_image[thread_1D_pos+2] = 255;
  106.         }
  107.         else if(i/64 == 1)
  108.         {
  109.                 d_image[thread_1D_pos] = 0;
  110.                 d_image[thread_1D_pos+1] = 255;
  111.                 d_image[thread_1D_pos+2] = 255-(i%64)*4;
  112.         }
  113.         else if(i/64 == 2)
  114.         {
  115.                 d_image[thread_1D_pos] = (i%64)*4;
  116.                 d_image[thread_1D_pos+1] = 255;
  117.                 d_image[thread_1D_pos+2] = 0;
  118.         } else {
  119.                 d_image[thread_1D_pos] = 255;
  120.                 d_image[thread_1D_pos+1] = 255-(i%64)*4;
  121.                 d_image[thread_1D_pos+2] = 0;
  122.         }
  123. }
  124.  
  125.  
  126. /*
  127. * Initialization of cuda
  128. */
  129. void operations(uchar* h_image, const int numRows, const int numCols, double cx, double cy, double zoom)
  130. {
  131.         int BLOCK_SIZE_X = 20;
  132.         int BLOCK_SIZE_Y = 20;
  133.     const dim3 blockSize(BLOCK_SIZE_X, BLOCK_SIZE_Y, 1);
  134.         const dim3 gridSize(numCols/BLOCK_SIZE_X+1, numRows/BLOCK_SIZE_Y+1, 1);
  135.  
  136.         uchar* d_image;
  137.         cudaMalloc(&d_image, sizeof(uchar) * 3 * numRows * numCols);
  138.  
  139.         // Menderlbrot
  140.         mandelbrot<<<gridSize, blockSize>>>(d_image, numRows, numCols, cx, cy, zoom, iterationMax);
  141.         cudaDeviceSynchronize();
  142.  
  143.         cudaMemcpy(h_image, d_image, sizeof(uchar) * 3 * numRows * numCols, cudaMemcpyDeviceToHost);
  144.         cudaFree(d_image);
  145. }
  146.  
  147.  
  148.  
  149.  
  150.  
  151.  
  152.  
  153. /*
  154. * Mouse callback
  155. */
  156. void mouseEvent(int evt, int x, int y, int flags, void* param){
  157.     if(evt==CV_EVENT_LBUTTONDOWN || evt == CV_EVENT_RBUTTONDOWN){
  158.  
  159.                 for(int i = 0 ; i < framesPerZoom ; i++) {
  160.                         // New center
  161.                         cx += (x-numCols/2)/zoom/framesPerZoom;
  162.                         cy += (y-numRows/2)/zoom/framesPerZoom;
  163.                         //printf("x=%d y=%d\ncx=%f, cy=%f zoom=%f\n", x, y, cx, cy, zoom);
  164.  
  165.                
  166.                         // New zoom
  167.                         zoom *= (evt==CV_EVENT_LBUTTONDOWN) ? std::pow(1.5, 1./framesPerZoom) : (evt==CV_EVENT_RBUTTONDOWN) ? std::pow(0.66, 1./framesPerZoom) : 1.0;
  168.  
  169.                         // New image
  170.                         //GpuTimer timer;
  171.                         //timer.Start();
  172.                         operations(h_image, numRows, numCols, cx, cy, zoom);
  173.                         //timer.Stop();
  174.                         //printf("Total time : %f ms\n", timer.Elapsed());
  175.  
  176.                         cv::imshow("Mandelbrot", cv::Mat(numRows, numCols, CV_8UC3, h_image));
  177.                         cv::waitKey(10);
  178.                 }
  179.     }
  180. }
  181.  
  182. /*
  183. * Main
  184. */
  185. int main(int argc, char **argv) {
  186.  
  187.        
  188.         h_image = (uchar*) malloc(sizeof(uchar) * 3 * numRows * numCols);
  189.  
  190.         GpuTimer timer;
  191.         timer.Start();
  192.         operations(h_image, numRows, numCols, cx, cy, zoom);
  193.         timer.Stop();
  194.         printf("Total time : %f ms\n", timer.Elapsed());
  195.  
  196.         cv::imshow("Mandelbrot", cv::Mat(numRows, numCols, CV_8UC3, h_image));
  197.         cvSetMouseCallback("Mandelbrot", mouseEvent, 0);
  198.  
  199.         cv::waitKey();
  200.         return 0;
  201. }
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
Top