CUDA 编程实例:计算点云法线

程序参考文章:http://blog.csdn.net/gamesdev/article/details/17535755  程序优化2

简介:CUDA ,MPI,Hadoop都是并行运算的工具。CUDA是基于NVIDIA GPU芯片计算。

阐述:GPU有很多个核(几百个),每个核可以跑一个线程,多个线程组成一个单位叫做块。


举个例子:
有三个向量 int a, b, c; 我们要计算a和b的向量之和存放到c中。

一般C语言:for(int i=0; i<10; i++)  c = a + b; 这个程序是顺序执行的!

CUDA编程做法:

GPU中的每个线程(核)有一个独立序号叫index,那么只要序号从0到9的线程执行c[index] = a[index] + b[index]; 就可以实现以上的for循环。

GPU的可贵之处就是,可以并发运行多个线程,相当于一个时间内赋值10次。

[cpp]   view plain  copy
  1.   
  2. cuda.cu  
  3.   
  4. #include <stdio.h>  
  5. #include <cuda_runtime.h>  
  6.   
  7. //运行在GPU  
  8. __global__ void vectorADD(int* a, int* b, int* c)  
  9. {  
  10.      int index = threadIdx.x; //获得当前线程的序号  
  11.      if(index < blockDim.x)  
  12.          c = a + b;  
  13. }  
  14. int main ()  
  15. {  
  16.         //定义10个GPU运算线程  
  17.         int N = 10;       
  18.         // 本地开辟三个数组存放我们要计算的内容  
  19.         int* h_a = (int*) malloc (N * sizeof(int));  
  20.         int* h_b = (int*) malloc (N * sizeof(int));  
  21.         int* h_c = (int*) malloc (N * sizeof(int));  
  22.         // 初始化数组A, B和C  
  23.         for(int i=0; i<N; i++) {  
  24.                 h_a = i;  
  25.                 h_b = i;  
  26.                 h_c = 0;  
  27.         }  
  28.         // 计算10个int型需要的空间  
  29.         int size = N * sizeof(int);     
  30.         // 在GPU上分配同样大小的三个数组  
  31.         int* d_a;  
  32.         int* d_b;  
  33.         int* d_c;  
  34.         cudaMalloc((void**)&d_a, size);  
  35.         cudaMalloc((void**)&d_b, size);  
  36.         cudaMalloc((void**)&d_c, size);  
  37.   
  38.         // 把本地的数组拷贝进GPU内存  
  39.         cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);  
  40.         cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);  
  41.         cudaMemcpy(d_c, h_c, size, cudaMemcpyHostToDevice);  
  42.           
  43.         // 定义一个GPU运算块 由 10个运算线程组成  
  44.         dim3 DimBlock = N;  
  45.         // 通知GPU用10个线程执行函数vectorADD  
  46.         vectorADD<<<1, DimBlock>>>(d_a, d_b, d_c);  
  47.         //  将GPU运算完的结果复制回本地  
  48.         cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);          
  49.         // 释放GPU的内存  
  50.         cudaFree(d_a);  
  51.         cudaFree(d_b);  
  52.         cudaFree(d_c);  
  53.         //  验证计算结果  
  54.         for(int j=0; j<N; j++)  
  55.                 printf("%d ", h_c);  
  56.         printf("\n");  
  57. }  

警告!:这个例子是编译不通过的;

首先:对 threadidx的使用,只能在CU文件里面;

其次:在cu文件里初始化数组是错误的: int * a ; a = new int [x];是错误的;  并且 malloc也是不可以的;

再者:文件路径里面不能包含中文,否则会出现 MSB8791 这种错误!


2. 利用CUDA并行计算点云法线

两个函数都存在于CU文件里! 通过外部CPP文件函数进行调用 

[cpp]   view plain  copy
  1. void normalEstimate(  
  2.     pcl::PointCloud<pcl::PointXYZRGB> &input ,  
  3.     pcl::PointCloud<pcl::PointXYZRGB> &output,  
  4.     int  k_,  
  5.     float search_parameter_,  
  6.     int THREAD_NUM  
  7.     )  

//运行在GPU//cal the Normal

[cpp]   view plain  copy
  1. __global__ void normalEstimateSingle(pcl::PointCloud<pcl::PointXYZRGB> &input ,pcl::PointCloud<pcl::PointXYZRGB> &output, int* nn_indices ,int*   nn_dists, int Gap, float search_parameter_)  
  2. {    
  3.     const size_t computeSize =input.size() / Gap;    
  4.     const size_t tID = size_t(threadIdx.x );           
  5.     int Mark;  
  6.     clock_t startTime; // 开始计时  
  7.     if ( tID == 0 ) startTime =clock( );// 选择任意一个线程进行计时       
  8.     //Thread loop!//循环发现邻域!寻找法线!  
  9.     for ( size_t idx = tID *computeSize; idx < ( tID + 1 ) * computeSize && idx < input.size(); ++idx )  {    
  10.       // pOut[threadIdx.x] += pIn[i] * pIn[i];    
  11.        Mark = pcl::searchForNeighbors (idx, search_parameter_, nn_indices, nn_dists);//对第IDX个建立索引!  
  12.        if (Mark == 0){  
  13.                 output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits<float>::quiet_NaN ();  
  14.                 continue;  
  15.             }  
  16.        else {  
  17.                 if (!isFinite (input[idx]) || Mark == 0){  
  18.                     output.points[idx].normal[0] = output.points[idx].normal[1] = output.points[idx].normal[2] = output.points[idx].curvature = std::numeric_limits<float>::quiet_NaN ();  
  19.                     continue;  
  20.                 }  
  21.             }  
  22.        pcl::computePointNormal (input, nn_indices,output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2], output.points[idx].curvature);  
  23.        pcl::flipNormalTowardsViewpoint (input_->points[idx], vpx_, vpy_, vpz_,  
  24.        output.points[idx].normal[0], output.points[idx].normal[1], output.points[idx].normal[2]);  
  25.   
  26.     }     
  27.     if ( tID == 0 ) *pElapsed =clock( ) - startTime;// 结束计时,返回至主程序    
  28. }  

//运行在CPU端!

[cpp]   view plain  copy
  1. // as the input  
  2. extern "C" void normalEstimate(  
  3.     pcl::PointCloud<pcl::PointXYZRGB> &input ,  
  4.     pcl::PointCloud<pcl::PointXYZRGB> &output,  
  5.     int  k_,  
  6.     float search_parameter_,  
  7.     int THREAD_NUM  
  8.     )  
  9. {  
  10.      // 在GPU上分配同样大小的三个数组  
  11.       pcl::PointCloud<pcl::PointXYZRGB> &inputX ;  
  12.       pcl::PointCloud<pcl::PointXYZRGB> &outputX;  
  13.        int* nn_indices ;  
  14.        int*   nn_dists;  
  15.          
  16.   
  17.     // 1、设置设备    
  18.     cudaError_t cudaStatus = cudaSetDevice( 0 );// 只要机器安装了英伟达显卡,那么会调用成功    
  19.     if ( cudaStatus != cudaSuccess )    
  20.     {    
  21.        fprintf( stderr, "调用cudaSetDevice()函数失败!" );    
  22.        return ;//false;    
  23.     }    
  24.   
  25.     // 使用CUDA内存分配器分配host端    
  26.     //cudaError_t cudaStatus = cudaMallocHost( &inputX, input.size() * sizeof( pcl::pointXYZRGB ) );    
  27.     //cudaError_t cudaStatus = cudaMallocHost( &outputX, output.size() * sizeof( pcl::Normal ) );    
  28.   
  29.     // 2、分配显存空间    
  30.     cudaError_t cudaStatus  = cudaMalloc( &inputX, input.size() * sizeof( pcl::pointXYZRGB ) );    
  31.     cudaError_t cudaStatusX = cudaMalloc( &outputX, output.size() * sizeof( pcl::Normal ) );  
  32.   
  33.       // cudaStatus = cudaMalloc( (void**)&pData, DataSize * sizeof( int) );    
  34.        if ( cudaStatus != cudaSuccess)    
  35.        {    
  36.            fprintf( stderr, "调用cudaMalloc()函数初始化显卡中数组时失败!" );    
  37.            break;    
  38.        }    
  39.   
  40.         // 3、将宿主程序数据复制到显存中    
  41.        cudaError_t cudaStatus2  = cudaMemcpy( inputX, input, input.size() * sizeof( pcl::pointXYZRGB ),cudaMemcpyHostToDevice );    
  42.        cudaError_t cudaStatusX2 = cudaMemcpy(outputX,output,output.size() * sizeof( pcl::pointXYZRGB ),cudaMemcpyHostToDevice );  
  43.        if ( cudaStatus != cudaSuccess)    
  44.        {    
  45.            fprintf( stderr, "调用cudaMemcpy()函数初始化宿主程序数据数组到显卡时失败!" );    
  46.            break;    
  47.        }    
  48.   
  49.        //cudaMalloc( (void**)&nn_dists,   k_ * sizeof( int) );    
  50.        //cudaMalloc( (void**)&nn_indices, k_ * sizeof( int) );    
  51.        //cudaMalloc( (void**)&Normal3f,  3 * sizeof( float) );     
  52.   
  53.      // 4、执行程序,宿主程序等待显卡执行完毕    
  54.        normalEstimateSingle<<<1, THREAD_NUM>>>( inputX,outputX, nn_indices, nn_dists, THREAD_NUM ,search_parameter_);  
  55.        //normalEstimateSingle(pcl::PointCloud<pcl::PointXYZRGB> &input ,pcl::PointCloud<pcl::PointXYZRGB> &output, int* nn_indices ,int*   nn_dists, int Gap)  
  56.          
  57.        // 5、查询内核初始化的时候是否出错    
  58.        cudaStatus = cudaGetLastError( );    
  59.        if ( cudaStatus != cudaSuccess)    
  60.        {    
  61.            fprintf( stderr, "显卡执行程序时失败!" );    
  62.            break;    
  63.        }  
  64.   
  65.      // 6、与内核同步等待执行完毕    
  66.        cudaStatus = cudaDeviceSynchronize( );    
  67.        if ( cudaStatus != cudaSuccess)    
  68.        {    
  69.            fprintf( stderr, "在与内核同步的过程中发生问题!" );    
  70.            break;    
  71.        }    
  72.      
  73.        // 7、获取数据  //只复制出法线即可!  
  74.        cudaStatus = cudaMemcpy(output,outputX,output.size() * sizeof( pcl::pointXYZRGB ),cudaMemcpyHostToDevice );    
  75.        if ( cudaStatus != cudaSuccess)    
  76.        {    
  77.            fprintf( stderr, "在将结果数据从显卡复制到宿主程序中失败!" );    
  78.            break;    
  79.        }    
  80.   
  81.        cudaFree( outputX );  
  82.        cudaFree(  inputX );  
  83. }  

注意事项:运行在GPU的函数,只能是原子函数,详情请见 《高性能并行编程实践》
  • 0
    点赞
  • 6
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值