《CUDA By Example》第四章的示例，觉得挺有趣的。由于里面绘图用到了cpu_bitmap.h头文件，又懒得去找书上的示例代码，就直接用opencv改了一下，直接将生成的Julia集存储到图像文件中。
1、Julia集
生成Julia集的算法十分简单。Julia集的基本算法是，通过一个简单的迭代等式对复平面中的点求值。如果在计算某个点时，迭代等式的结果是分散的，那么这个点就不属于Julia集合。相反，如果迭代等式计算得到的一系列值都位于某个边界范围内，那么这个点就属于Julia集合。迭代等式如下：
![][1]
[1]: http://latex.codecogs.com/svg.latex?Z_{n+1}=Z_n^2+C
迭代过程包括计算当前值的平方，然后在加上一个常数C得到下一个值。
2、复数结构体cuComplex的定义
结构体包含一个构造函数，一个计算复数模的平方函数magnitude2()以及根据复数运算规则对乘法和加法的运算符重载函数。
![][2]
[2]: http://latex.codecogs.com/svg.latex?(a+bi)*(c+di)=(ac-bd)+(ad+bc)i
![][3]
[3]: http://latex.codecogs.com/svg.latex?(a+bi)+(c+di)=(a+c)+(b+d)i
在《CUDA By Example》书中，构造函数没有添加__device__前缀，程序无法通过编译。但是后来也发现了别的问题，当图片边长过大，程序运行不会报错但是无法正确生成Julia集图片，怀疑是因为在GPU上迭代计算时生成了过多cuComplex对象导致内存不足。但是也只是猜测，自己对C++内存回收机制也不是了解，等有空了仔细学学。
struct cuComplex{
float r;
float i;
__device__ 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);
}
};
3、Julia集判断函数
该函数功能为判断坐标(i,j)是否属于Julia集合。首先将像素坐标转换为复数空间的坐标，(jx, jy)则为转换后的复数空间坐标，计算出复数空间坐标后，迭代计算判断是否属于Julia集合，其中c是一个复数常量，当选择不同的值时，可以生成不同的图片。
该函数中计算了100次迭代，每次迭代完成后都判断结果是否超过阈值1000，如果超过则说明不属于Julia集，返回false；如果100次迭代完成后都没有返回false则返回true说明该坐标属于Julia集。
__device__ bool 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 c(0.285, 0.02);
cuComplex a(jx, jy);
for(int i = 0; i < 100; i++)
{
a = a*a+c;
if(a.magnitude2() > 1000)
return false;
}
return true;
}
4、核函数和main函数
在main函数中制定了多个线程块来执行函数kernel()，申明了一个二维的线程格grid
dim3 grid(DIM, DIM);
其中DIM为生成的图片边长，也就是长和宽上的像素数；然后将dim3变量grid传递给CUDA运行:
kernel<<<grid, 1>>>(d_img, DIM);
d_img为设备上分配的内存空间，调用kernel函数将对其进行修改，将其中属于Julia集合的像素坐标设置为绿色，否则设置为黑色。
__global__ void kernel(uchar4 *d_img, int dim)
{
//获取一维索引
int i = blockIdx.x;
int j = blockIdx.y;
int offset = i+j*dim;
if(julia(i, j)){
d_img[offset].x = 0;
d_img[offset].y = 255;
d_img[offset].z = 0;
}
else
{ d_img[offset].x = 0;
d_img[offset].y = 0;
d_img[offset].z = 0;
}
}
int main()
{
Mat img(DIM, DIM, CV_8UC3);
uchar4 *d_img, *i_img;
i_img = (uchar4*)malloc(DIM*DIM*sizeof(uchar4));
cudaMalloc(&d_img, DIM*DIM*sizeof(uchar4));
dim3 grid(DIM, DIM);
kernel<<<grid, 1>>>(d_img, DIM);
cudaMemcpy(i_img, d_img, DIM*DIM*sizeof(uchar4),cudaMemcpyDeviceToHost);
for(int i = 0; i < DIM; ++i)
{
for(int j = 0; j < DIM; ++j)
{
img.at<Vec3b>(i, j)[0] = i_img[i*DIM+j].x;
img.at<Vec3b>(i, j)[1] = i_img[i*DIM+j].y;
img.at<Vec3b>(i, j)[2] = i_img[i*DIM+j].z;
}
}
imwrite("img.jpg", img);
//内存释放
cudaFree(d_img);
free(i_img);
return 0;
}
5、完整代码如下：
//kernel.cu
//author:Curya
//Date:2017-07-03
#include <iostream>
#include <cuda_runtime.h>
#include <opencv2\core\core.hpp>
#include <opencv2\imgproc\imgproc.hpp>
#include <opencv2\highgui\highgui.hpp>
#define DIM 600
using namespace cv;
struct cuComplex{
float r;
float i;
__device__ 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);
}
};
__device__ bool 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 c(0.285, 0.02);
cuComplex c(0, 0.73);
cuComplex a(jx, jy);
for(int i = 0; i < 100; i++)
{
a = a*a+c;
if(a.magnitude2() > 1000)
return false;
}
return true;
}
__global__ void kernel(uchar4 *d_img, int dim)
{
int i = blockIdx.x;
int j = blockIdx.y;
int offset = i+j*dim;
if(julia(i, j)){
d_img[offset].x = 0;
d_img[offset].y = 255;
d_img[offset].z = 0;
}
else
{ d_img[offset].x = 0;
d_img[offset].y = 0;
d_img[offset].z = 0;
}
}
int main()
{
Mat img(DIM, DIM, CV_8UC3);
uchar4 *d_img, *i_img;
i_img = (uchar4*)malloc(DIM*DIM*sizeof(uchar4));
cudaMalloc(&d_img, DIM*DIM*sizeof(uchar4));
dim3 grid(DIM, DIM);
kernel<<<grid, 1>>>(d_img, DIM);
cudaMemcpy(i_img, d_img, DIM*DIM*sizeof(uchar4),cudaMemcpyDeviceToHost);
for(int i = 0; i < DIM; ++i)
{
for(int j = 0; j < DIM; ++j)
{
img.at<Vec3b>(i, j)[0] = i_img[i*DIM+j].x;
img.at<Vec3b>(i, j)[1] = i_img[i*DIM+j].y;
img.at<Vec3b>(i, j)[2] = i_img[i*DIM+j].z;
}
}
imwrite("img.jpg", img);
//内存释放
cudaFree(d_img);
free(i_img);
return 0;
}
