CUDA实战笔记(5)——绘制Julia集

本系列代码托管于：https://github.com/chintsan-code/cuda_by_example
本篇使用的项目为：julia_cpu、julia_gpu

上一篇笔记实现了GPU上的矢量相加，本篇实现绘制Julia集。对于Julia集，只需要理解成满足某个公式的得所有点构成的集合即可，不用深入理解，这个公式是：

$$Z_{n+1} = Z_n^2 + C$$

可以看到，这个公式的迭代过程为：首先计算当前值得平方，然后再加上一个常数，就得到等式的下一个值。

一. 基于CPU的Julia集绘制

// julia_cpu

#include <stdio.h>
#include "../../common/cpu_bitmap.h"

#define DIM 1000

struct cuComplex
{
    float r;  // 复数的实数部分
    float i;  // 复数的虚数部分

    cuComplex(float a, float b) :r(a), i(b) { }

    float magnitude2(void)
    {
        return r * r + i * i;  // 复数的模的平方
    }

    cuComplex operator * (const cuComplex& a)
    {
        return cuComplex(r * a.r - i * a.i, i * a.r + r * a.i);
    }

    cuComplex operator + (const cuComplex& a)
    {
        return cuComplex(r + a.r, i + a.i);
    }
};

int julia(int x, int y)
{
    const float scale = 1.5;
    // DIM / 2 - x、DIM / 2 - y将原点定位到图像中心
    // 除以(DIM / 2)是为了确保图像的范围为[-1.0, 1.0]
    // scale是用来缩放图像的，可以自行修改
    float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);
    float jy = scale * (float)(DIM / 2 - y) / (DIM / 2);

    cuComplex c(-0.8, 0.156);  // 复数C = -0.8 + 0.156i
    cuComplex z(jx, jy);

    int i = 0;
    for (i = 0; i < 200; i++)
    {
        z = z * z + c; //Zn+1 = Zn^2 + C
        if (z.magnitude2() > 1000)
        {
            // 迭代200次，每次迭代完都判断结果是否超过阈值(这里是1000)，如果超过就不属于julia集
            return 0;
        }
    }
    return 1;  // 属于Julia集
}

void kernel(unsigned char* ptr)
{
    for (int y = 0; y < DIM; y++)
    {
        for (int x = 0; x < DIM; x++)
        {
            int offset = x + y * DIM;  // 像素在内存中的线性偏移，因为图像在内存中实际是一维存储的

            int juliaValue = julia(x, y);  // 判断点(x, y)是否属于Julia集合，属于返回1，不属于返回0
            // juliaValue为0时为黑色(0,0,0)，为1时为红色(255,0,0)
            ptr[offset * 4 + 0] = 255 * juliaValue; // red通道
            ptr[offset * 4 + 1] = 0;                // green通道
            ptr[offset * 4 + 2] = 0;                // blue通道
            ptr[offset * 4 + 3] = 255;              // alpha通道
        }
    }
}

int main()
{
    CPUBitmap bitmap(DIM, DIM);
    unsigned char* ptr = bitmap.get_ptr();
    kernel(ptr);  // 将指向图像的指针传递给核函数
    bitmap.display_and_exit();

    getchar();
    return 0;
}

首先看主函数main：

int main()
{
    CPUBitmap bitmap(DIM, DIM);
    unsigned char* ptr = bitmap.get_ptr();
    kernel(ptr);  // 将指向图像的指针传递给核函数
    bitmap.display_and_exit();

    getchar();
    return 0;
}

通过工具库创建一张位图，大小是DIM * DIM。然后获取该图像的指针，将其传递给核函数，最后将图像显示处理。

然后是核函数kernel:

void kernel(unsigned char* ptr)
{
    for (int y = 0; y < DIM; y++)
    {
        for (int x = 0; x < DIM; x++)
        {
            int offset = x + y * DIM;

            int juliaValue = julia(x, y);  // 判断点(x, y)是否属于Julia集合，属于返回1，不属于返回0
            // juliaValue为0时为黑色(0,0,0)，为1时为红色(255,0,0)
            ptr[offset * 4 + 0] = 255 * juliaValue; // red通道
            ptr[offset * 4 + 1] = 0;                // green通道
            ptr[offset * 4 + 2] = 0;                // blue通道
            ptr[offset * 4 + 3] = 255;              // alpha通道
        }
    }
}

核函数对图像中的所有点进行后台绘制，一共有DIM * DIM个点，其中会调用julia函数对当前点进行判断，判断其是否属于Julia集，如果属于，将这个点填充为红色；如果不属于，填充为黑色。

julia函数和cuComplex结构体都是为了Julia集服务了，书中有介绍，这里不写了，仅对关键部分做了注释。

二. 基于GPU的Julia集绘制

// julia_gpu

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include "../../common/book.h"
#include "../../common/cpu_bitmap.h"

#define DIM 1000

struct cuComplex
{
    float   r;
    float   i;
    // cuComplex( float a, float b ) : r(a), i(b)  {}
    __device__ cuComplex(float a, float b) : r(a), i(b) { }  // __device__表示代码将在设备而不是主机上运行。
                                                             // 由于声明为__device__函数，因此只能从其他__device__
                                                             // 或__global__函数中调用他们
    __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__ 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;
}

__global__ void kernel(unsigned char* ptr) {
    // 将threadIdx/blockIdx映射到像素位置
    int x = blockIdx.x;
    int y = blockIdx.y;

    int offset = x + y * gridDim.x;  // 对所有的线程块而言，gridDim是常数，用来保存线程格每一维的大小，
                                     // gridDim.x是线程格的宽度

    int juliaValue = julia(x, y);
    ptr[offset * 4 + 0] = 255 * juliaValue;
    ptr[offset * 4 + 1] = 0;
    ptr[offset * 4 + 2] = 0;
    ptr[offset * 4 + 3] = 255;
}

int main() {
    CPUBitmap bitmap(DIM, DIM);
    unsigned char* dev_ptr;

    HANDLE_ERROR(cudaMalloc((void**)&dev_ptr, bitmap.image_size()));

    dim3 grid(DIM, DIM);  // 二维线程格(Grid)。dim3表示一个三维数组，在这里第三维其实就是1，同dim3 grid(DIM, DIM, 1)
    kernel<<<grid, 1>>>(dev_ptr);

    HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(), dev_ptr, bitmap.image_size(), cudaMemcpyDeviceToHost));

    bitmap.display_and_exit();

    HANDLE_ERROR(cudaFree(dev_ptr));

    return 0;
}

由于在这个绘制情景下，由于图像中的每个点的计算都是独立的，计算结果与其他点无关，因此可以指定多个线程并行执行kernel函数。而且在这个图像处理的问题中，Block的二维索引会带来非常大的便利。

首先声明了一个三维的线程格(Grid)：

dim3 grid(DIM, DIM);

类型dim3并不是标准C中定义的类型，它被定义在CUDA相关的头文件中，它表示一个三维数组，可以用于指定启动的Block的数量，在这里第三维其实就是1。因此下面的代码表示告诉CUDA Runtime，我们需要创建一个二维的Grid，其中包含DIM * DIM个Block。

dim3 grid(DIM, DIM);
kernel<<<grid, 1>>>(dev_ptr);

回到kernel函数：

__global__ void kernel(unsigned char* ptr) {
    int x = blockIdx.x;
    int y = blockIdx.y;

    int offset = x + y * gridDim.x;  // 对所有的线程块而言，gridDim是常数，用来保存线程格每一维的大小，
                                     // gridDim.x是线程格的宽度

    int juliaValue = julia(x, y);
    ptr[offset * 4 + 0] = 255 * juliaValue;
    ptr[offset * 4 + 1] = 0;
    ptr[offset * 4 + 2] = 0;
    ptr[offset * 4 + 3] = 255;
}

与CPU版本不同的是，不再需要双重for循环来逐个生成图像的每一个像素的索引，CUDA Runtime将在变量blockIdx中包含这些像素索引，由于我们设置了合适的Grid和Block，因此每一个像素都可以获得一个Block进行渲染绘制。

接下来需要得到ptr中的线性偏移(与基于CPU的一致，因为图像在内存中实际是一维存储的)。这个偏移值通过内置变量gridDim来得到，实际上，但我们告诉CUDA Runtime需要创建一个包含DIM*DIM个Block的二维Grid时，这个 gridDim变量就已经被确定了，相当于常数，大小为(DIM, DIM)

julia函数和cuComplex结构体和CPU版本的基本一致，二者的差异主要在多了修饰符__device__，这表示该段代码将在GPU而不是主机上运行。由于这些函数已声明为__device__函数，因此只能从其他的__device__函数或__global__函数中调用它们。

参考：

• 《GPU高性能编程 CUDA实战》