詳解 CUDA By Example 中的 Julia Set 繪製GPU優化

筆者測試環境VS2019。

基本介紹

原書做者引入Julia Sets意在使用GPU加速圖形的繪製。Julia Set 是指知足下式迭代收斂的複數集合
\[ Z_{n+1}=Z_{n}^2+C \]

環境配置

跑這個例子的主要困難應該在於配置環境。這個程序依賴於openGL中的glut庫。因爲VS2019的整個軟件架構發生了很大變化，一些連接庫和頭文件的位置都發生了改變，所以一些文章中的配置方法失效了。

首先咱們須要獲取glut庫的頭文件以及動態連接庫。

點擊這裏cg-toolkit獲取。安裝成功以後，找到C:\Program Files (x86)\NVIDIA Corporation\Cg。注意勾選安裝選項的x64相關應用。

將其中的lib文件夾中的_glut32.lib_複製到C:\Program Files (x86)\Windows Kits\10\Lib\10.0.18362.0\ucrt\x86

將其中的lib.x64文件夾中的glut32.lib複製到C:\Program Files (x86)\Windows Kits\10\Lib\10.0.18362.0\ucrt\x64而且重命名其爲glut64.lib

筆者運行的是64位系統，就將bin.x64中的_glut32.dll_複製到C:\Windows\System32下

在這裏下載頭文件。下載完成以後，將頭文件拷貝到C:\Program Files (x86)\Windows Kits\10\Include\10.0.18362.0\ucrt。並創建文件夾GL把它們包括起來。

提示，核心是找到C:\Program Files (x86)\Windows Kits\10，不要在Microsoft Visual Studio文件夾裏浪費時間。

後面的10.0.18362.0根據版本不一樣可能不一致，具體問題具體分析

這個代碼還須要一些別的頭文件。如gl_helper.h, book.h, cpu_bitmap.h 等 在這裏下載後複製到C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v10.2\include

主要代碼

CPU Julia Set

RGBA模式中，每個像素會保存如下數據：R值（紅色份量）、G值（綠色份量）、B值（藍色份量）和A值（alpha份量）。其中紅、綠、藍三種顏色相組合，就能夠獲得咱們所須要的各類顏色，而alpha不直接影響顏色，它的含義是透明度。1

下面是純粹CPU中的代碼，基本的註釋在代碼中

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "device_functions.h"
#include "device_atomic_functions.h"
#include <cuda.h>
#include "book.h"
#include <cpu_bitmap.h>
#include <stdio.h>

#define DIM 1000        //圖像的像素邊長大小

struct cuComplex
{
    float r;
    float i;
    cuComplex(float a, float b) : r(a), i(b) {}
    float magnitude2() { 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;    //放大倍率
    float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);    //座標變換，投影到-1~1scale
    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;
}

void kernel(unsigned char* ptr)
{
    for (int y = 0; y < DIM; y++)   //遍歷整個bitmap
    {
        for (int x = 0; x < DIM; x++)
        {
            int offset = x + y * DIM;
            int juliaValue = julia(x, y);
            //注意openGL這裏的顏色格式是RGBA,000爲黑色
            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* ptr = bitmap.get_ptr();
    kernel(ptr);    //運行渲染
    bitmap.display_and_exit();
}

GPU Julia Set

注意因爲內核函數是global的，要在GPU上運行須要將其調用的julia函數加上device。又由於，device函數只能由device函數或者global函數調用，因此最好把結構體中的全部函數都加上device。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "device_functions.h"
#include "device_atomic_functions.h"
#include <cuda.h>
#include "book.h"
#include <cpu_bitmap.h>
#include <stdio.h>

//小於65536
#define DIM 1000        //圖像的像素邊長大小

struct cuComplex
{
    float r;
    float i;
    __device__ cuComplex(float a, float b) : r(a), i(b) {}
    __device__ float magnitude2() { 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);    //座標變換，投影到-1~1scale
    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)
{
    int x = blockIdx.x; //縱向線程索引(x方向朝右，是行)
    int y = blockIdx.y; //縱向線程索引(y方向朝下，是列）
    int offset = x + y * 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_bitmap;
    //在GPU中分配空間
    HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));
    dim3 grid(DIM, DIM);    //dim3結構體
    kernel <<<grid, 1 >>> (dev_bitmap); //一個線程塊中的線程網絡1000x1000
    HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(), dev_bitmap, bitmap.image_size(), cudaMemcpyDeviceToHost));    //將dev_bitmap中的內容從device拷貝到cpu中
    bitmap.display_and_exit();
    HANDLE_ERROR(cudaFree(dev_bitmap));
}

---

參考資料