Feb 212011
 

My “hello world” for a new language is to generate fractals and visualise the result. Having done this generally means I’m starting to grasp some of the principles.

JuliaSet_Result_small

In C++, I like to use the GraphicsMagick libraries to read and generate images. They are not perfect and tend not to scale very well (I would use DevIL in the scaling case) but they are easy to use. They enable me to paint a canvas and easily output the result to a .png (or any other) image.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
// Includes
#include <stdio.h>
#include <cuda.h>
#include <cuda_runtime_api.h>
#include <Magick++.h>
 
// What size fractal do I want?
#define DIM 1000
 
using namespace std;
using namespace Magick;
 
 
/*
 *  A struct: CUDA style
 *  A method for storing a complex number with single-precision
 *  floating point components
 */	   
struct cuComplex
{
   float r;
   float i;
 
   cuComplex( float a, float b) : r(a), i(b) {}
   __device__ float magnitude2( void )
   {
      return r*r+i*i;
   }
 
   __device__ cuComplex operator*(const cuComplex& a)
   {
      return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
   }
 
   __device__ cuComplex operator+(const cuComplex& a)
   {
      return cuComplex(r+a.r, i+a.i);
   }
};
 
 
/*
 *  Save the image
 */
void SaveImage(unsigned char *ptr)
{
   fprintf(stdout, "Writing Result Image (result.png)n");
   ColorGray c("white");
   Image image(Geometry(DIM, DIM), "black");
   image.type(GrayscaleType);
   for (int j = 0; j<DIM; j++)
   {
      for (int i = 0; i<DIM; i++)
      {		
         if (ptr[j*DIM + i] != 0)
            image.pixelColor(i, j, c );
      }
   }
   image.write("result.png");
}
 
 
/*
 *  return 1 if point inside the set, otherwise, 0
 */
__device__ int julia(int x, int y)
{
   const float scale = 1.5;
   float jx = scale * (float)(DIM/2 - x)/(DIM/2);
   float jy = scale * (float)(DIM/2 - y)/(DIM/2);
 
   cuComplex c(-0.8, 0.156);
   cuComplex a(jx, jy);
 
   int i = 0;
   for (i = 0; i < 200; i++)
   {
      a = a*a + c;
      if (a.magnitude2() > 1000)
         return 0;
   }
   return 1;
}
 
 
/*
 * A few new things here, see explaination below
 */
__global__ void kernel(unsigned char *ptr)
{
   // map from threadIdx/BlockIdx to pixel position
   int x = blockIdx.x;
   int y = blockIdx.y;
   int offset = x + y*gridDim.x;
 
   // now calculate the value at that position
   int juliaValue = julia(x,y);
   ptr[offset] = juliaValue;
}
 
 
/*
 *  The same-same Main
 */
int main( void)
{
   unsigned char *canvas;
   unsigned char *dev_canvas;
 
   //Allocate CPU memory
   canvas = new unsigned char[DIM*DIM];
 
   //Allocate GPU memory
   cudaMalloc((void**)&dev_canvas, DIM*DIM*sizeof(unsigned char));
 
   // See note below
   dim3 grid(DIM,DIM);
   kernel<<<grid,1>>>(dev_canvas);
 
   cudaMemcpy(canvas, dev_canvas, DIM*DIM*sizeof(unsigned char), cudaMemcpyDeviceToHost);
 
   cudaFree(dev_canvas);
 
   SaveImage(canvas);
 
   return 0;  
}

 Leave a Reply

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>

(required)

(required)

Human Conf Test * Time limit is exhausted. Please reload CAPTCHA.