利用CUDA和opencv绘制Julia集

《CUDA By Example》第四章的示例，觉得挺有趣的。由于里面绘图用到了cpu_bitmap.h头文件，又懒得去找书上的示例代码，就直接用opencv改了一下，直接将生成的Julia集存储到图像文件中。

---

1、Julia集

生成Julia集的算法十分简单。Julia集的基本算法是，通过一个简单的迭代等式对复平面中的点求值。如果在计算某个点时，迭代等式的结果是分散的，那么这个点就不属于Julia集合。相反，如果迭代等式计算得到的一系列值都位于某个边界范围内，那么这个点就属于Julia集合。迭代等式如下：

$$Z_{n+1}=Z_n^2+C$$
迭代过程包括计算当前值的平方，然后在加上一个常数C得到下一个值。

2、复数结构体cuComplex的定义

结构体包含一个构造函数，一个计算复数模的平方函数magnitude2()以及根据复数运算规则对乘法和加法的运算符重载函数。

$$(a+bi)*(c+di)=(ac-bd)+(ad+bc)i$$
$$(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;
}

6、运行结果

$C=-0.8+0.156i$

$C=0.285+0.02i$

$C=0.285+0.01i$

$C=-0.70176-0.3842i 迭代次数：25$