基于前两篇博客，其实我们对gpu编程已经掌握得差不多了，在这第三篇博客中，最要是两个例子，一个是光线追踪，一个是热传导的模拟。进一步介绍两中内存，constant memory和 texture memory。
光线追踪
光线追踪最是求三维场景在二维平面的投影，追踪三维物体的光线在平面的位置。这里我们介绍最简单的方法。
从成像平面的每一个像素出发，追踪它最终击中的三维物体的，然后根据击中点来绘制像素，这样遍历整个成像平面后，我们就得到了三维空间的成像图。原理如下图所示：
ray_tracer_principle.png
程序如下：
#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;
}
程序说明：
- 成像平面是x-y平面。
- 我们先声明了一个球类，这个类包含了球的空间坐标和半径，以及其颜色，另外我们定义了一个击中函数，来测试某个像素点发出的“逆光线”是否击中了球体。n是击中点离球心平面的一个度量，n越大，我们击中的位置越靠近球的投影中心，主要用来确定成像片面的像素的颜色明暗程度。返回值是击中点与成像平面的距离的第一度量。越大的话就越靠近成像平面，当一条“逆光线”击中不只一个球体的时候，我们应该选择最前面的点（最近的那个点）来绘制。未击中任何球，返回负无穷大。
- kernel函数中，if( t>maxz) 用来判定每次绘制都是最前的像素。其他的地方与之前的程序中大同小异。
constant Sphere s[SPHERE];声明一个常显存。
cudaMemcpyToSymbol(s, temp_s, SPHERE *sizeof(Sphere));是将内存中的数据拷贝至显存中。然后是计算。
本程序中多了cudaEvent_t这个变量，从名字可以看出，这是一个cuda的事件类型，我们这里主要用来测算cuda的运行性能。使用也很简单。
最终，在笔者的电脑是，采用constant类型的程序是不采用的程序耗时的7/10。这里不明显时候因为这里的constant类数据不大，当数据较大时，长内存的优势就会明显得多。
raay_tracer.png
热传递的模拟
热传递的模拟过成，我们在一幅图上，每个像素点只考虑它上下左右四个位置的强度，然后采用公式：
必须说明的是，这不是一个准确的热传递公式，甚至算不上是近似公式。其中的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的参数。关于这一点，笔者之前没注意到，随便写的参数，导致跑出的程序结果不符合预期。下面是运行的结果：
heat_emulator.png
网友评论