GPU编程三瞥

基于前两篇博客,其实我们对gpu编程已经掌握得差不多了,在这第三篇博客中,最要是两个例子,一个是光线追踪,一个是热传导的模拟。进一步介绍两中内存,constant memory和 texture memory。

光线追踪

光线追踪最是求三维场景在二维平面的投影,追踪三维物体的光线在平面的位置。这里我们介绍最简单的方法。

从成像平面的每一个像素出发,追踪它最终击中的三维物体的,然后根据击中点来绘制像素,这样遍历整个成像平面后,我们就得到了三维空间的成像图。原理如下图所示:

程序如下:

#include<opencv2/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<memory>
#include "cuda.h"
#include<iostream>
#include<memory>

using namespace std;
#define DIM 1024
#define INF 2e10f
#define SPHERE 20
#define rnd(x) (x*rand()/RAND_MAX)

class Sphere{
public:
    double r, g, b;
    double x, y, z, radius;
    __device__ float hit(double ox, double oy, double *n){
        double dx = ox-x;
        double dy = oy-y;
        double xy_2 = dx*dx + dy*dy;
        double rad_2 = radius*radius;
        if( xy_2 < radius*radius){
            float dz = sqrtf(rad_2 - xy_2);
            *n = dz/radius;
            return dz+z;
        }
        return -INF;
    }
};

__constant__ Sphere s[SPHERE];
__global__ void kernel(unsigned char *ptrs){
    int x =threadIdx.x + blockDim.x * blockIdx.x;
    int y =threadIdx.y + blockDim.y * blockIdx.y;
    int offset = x + y*blockDim.x * gridDim.x;

    double ox= x-DIM/2.0;
    double oy= y-DIM/2.0;
    double r=0, g=0, b=0; 
    double maxz = -INF;
    for(int i=0; i<SPHERE; i++){
        double n;
        double t = s[i].hit(ox, oy, &n);
        //if current hit is more closer to camera than last hit, use current data
        if( t>maxz){
            float fscale = n;
            r = s[i].r * fscale;
            g = s[i].g * fscale;
            b = s[i].b * fscale;
            maxz = t;//
        }
    }
    ptrs[4 * offset + 0] = int(r*255);
    ptrs[4 * offset + 1] = int(g*255);
    ptrs[4 * offset + 2] = int(b*255);
    ptrs[4 * offset + 3] = 255;
}

int main(){
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start, 0);
    cv::Mat_<cv::Vec3b> img(DIM,DIM);
    unsigned char ptrs[4*DIM*DIM];
    unsigned char *dev_ptrs;
    Sphere temp_s[SPHERE];

    srand((unsigned) time(0));
    for(int i=0; i<SPHERE; i++){
        temp_s[i].r = rnd(1.0f);
        temp_s[i].g = rnd(1.0f);
        temp_s[i].b = rnd(1.0f);
        temp_s[i].x = rnd(1000.0f)-500;
        temp_s[i].y = rnd(1000.0f)-500;
        temp_s[i].z = rnd(1000.0f)-500;
        temp_s[i].radius = rnd(100.0f) + 20;
    }
    cudaMalloc((void**)&dev_ptrs, 4*DIM*DIM*sizeof(unsigned char));
    //cudaMalloc((void**)&s, SPHERE*sizeof(Sphere));
    //cudaMemcpy(s, temp_s, SPHERE*sizeof(Sphere), cudaMemcpyHostToDevice);
    cudaMemcpyToSymbol(s, temp_s, SPHERE * sizeof(Sphere));

    dim3 blocks(32, 32);
    dim3 threads(32, 32);
    kernel<<<blocks, threads>>>(dev_ptrs);
    cudaMemcpy(ptrs, dev_ptrs, 4*DIM*DIM*sizeof(unsigned char), cudaMemcpyDeviceToHost);

    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    float mytime;
    cudaEventElapsedTime(&mytime, start, stop);
    cout<<"performace:\n"<<mytime<<endl;

    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    cudaFree(dev_ptrs);
    cudaFree(s);
    for(int y=0; y<img.rows; y++){
        for(int x=0; x<img.cols; x++){
            for(int ch=0; ch<3; ch++){
                img.at<cv::Vec3b>(x, y)[ch] = ptrs[ 4*(x+y*DIM) + ch];
            }
        }
    }
    cv::imshow("show", img);
    cv::waitKey(0);
    return 0;
}

程序说明:

  1. 成像平面是x-y平面。
  2. 我们先声明了一个球类,这个类包含了球的空间坐标和半径,以及其颜色,另外我们定义了一个击中函数,来测试某个像素点发出的“逆光线”是否击中了球体。*n是击中点离球心平面的一个度量,*n越大,我们击中的位置越靠近球的投影中心,主要用来确定成像片面的像素的颜色明暗程度。返回值是击中点与成像平面的距离的第一度量。越大的话就越靠近成像平面,当一条“逆光线”击中不只一个球体的时候,我们应该选择最前面的点(最近的那个点)来绘制。未击中任何球,返回负无穷大。
  3. kernel函数中,if( t>maxz) 用来判定每次绘制都是最前的像素。其他的地方与之前的程序中大同小异。
  4. __constant__ Sphere s[SPHERE];声明一个常显存。

    cudaMemcpyToSymbol(s, temp_s, SPHERE *sizeof(Sphere));是将内存中的数据拷贝至显存中。然后是计算。
  5. 本程序中多了cudaEvent_t这个变量,从名字可以看出,这是一个cuda的事件类型,我们这里主要用来测算cuda的运行性能。使用也很简单。

最终,在笔者的电脑是,采用constant类型的程序是不采用的程序耗时的7/10。这里不明显时候因为这里的constant类数据不大,当数据较大时,长内存的优势就会明显得多。

热传递的模拟

热传递的模拟过成,我们在一幅图上,每个像素点只考虑它上下左右四个位置的强度,然后采用公式:

$$C_{n+1}=T_n+\sum_{neighbor} k \dot (T_{neighbor}-T_n)$$

必须说明的是,这不是一个准确的热传递公式,甚至算不上是近似公式。其中的k表征传递的速度。程序如下:

#include<opencv2/core.hpp>
#include<opencv2/highgui/highgui.hpp>
#include<memory>
#include "cuda.h"
#include<iostream>
#include<memory>

using namespace std;
#define DIM 1024
//系数应当小于0.25
#define SPEED 0.25

__global__ void mycopy(const float *src, float *dest){
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;
    if( src[offset]!=0) dest[offset] = src[offset];
    //dest[offset] = src[offset];
}

__global__ void kernel(float *optrs, const float *iptrs, int tick){
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    int right = offset + 1;
    int left = offset - 1;
    if(x==DIM-1) right--;
    if(x==0) left++;

    int top = offset - DIM;
    int down = offset + DIM;
    if(y==DIM-1) down -= DIM;
    if(y==0) top += DIM;

    optrs[offset] = iptrs[offset] + SPEED*(iptrs[right] + iptrs[left] + iptrs[top] + iptrs[down] - 4*iptrs[offset] );
}

int main(){
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start, 0);

    cv::Mat_<cv::Vec3b> img(DIM,DIM);
    float img_ptrs[DIM*DIM]={0};
    float *out_ptr;
    float *const_ptr;
    float *in_ptr;
    dim3 blocks(32, 32);
    dim3 threads(32, 32);

    cudaMalloc((void**)&out_ptr, DIM*DIM*sizeof(float));
    cudaMalloc((void**)&const_ptr, DIM*DIM*sizeof(float));
    cudaMalloc((void**)&in_ptr, DIM*DIM*sizeof(float));
    
    for(int i=0; i<DIM; i++){
        for(int j=0; j<DIM; j++){
            if(i>110 && i<210 && j>110 && j<210)
                img_ptrs[i*DIM + j] = 255;
            if(i>210 && i<310 && j>210 && j<310)
                img_ptrs[i*DIM + j] = 255;
        }
    }

    cudaMemcpy( const_ptr, img_ptrs, DIM*DIM*sizeof(float), cudaMemcpyHostToDevice );

    for(int i=0; i<DIM*DIM; i++)
        img_ptrs[i]=70;
    cudaMemcpy( in_ptr, img_ptrs, DIM*DIM*sizeof(float), cudaMemcpyHostToDevice );

    for(int it=0; it<9000; it++){
        mycopy<<<blocks,threads>>>(const_ptr, in_ptr);

        kernel<<<blocks, threads>>>(out_ptr, in_ptr, it);
        swap(in_ptr, out_ptr);
        //cudaMemcpy(img_ptrs, in_ptr, DIM*DIM*sizeof(float), cudaMemcpyDeviceToHost);
        cudaMemcpy(img_ptrs, in_ptr, DIM*DIM*sizeof(float), cudaMemcpyDeviceToHost);

        cout<<it<<" pixel(215,205): "<<img_ptrs[215 + 205*DIM]<<endl;
        for(int i=0; i< img.rows; i++){
            for(int j=0; j<img.cols; j++){
                for(int ch=0; ch<3; ch++)
                    img.at<cv::Vec3b>(i,j)[ch]=img_ptrs[ j*DIM+i];
            }
        }
        if(it == 8999)
            cv::waitKey(0);
        cv::imshow("test", img);
        cv::waitKey(1);
    }

    cudaEventRecord(stop, 0);
    cudaEventSynchronize(stop);
    float mytime;
    cudaEventElapsedTime(&mytime, start, stop);
    cout<<"performace:\n"<<mytime<<endl;

    cudaEventDestroy(start);
    cudaEventDestroy(stop);
    cudaFree(in_ptr);
    cudaFree(out_ptr);
    cudaFree(const_ptr);

    return 0;
}

程序说明:

特别强调一点,程序中SPEED应当是大于零,小于0.25的参数。关于这一点,笔者之前没注意到,随便写的参数,导致跑出的程序结果不符合预期。下面是运行的结果:

 

发表评论

电子邮件地址不会被公开。 必填项已用*标注