用cuda实现图像缩放(从中领悟到了学习cuda编程的方法)

         最近在cuda实现HOG特征抽取。感觉算法中有不少地方可以并行化，但是怎么并行化才会优化性能，自己还没有找到很明确的方法。
       HOG特征抽取有一个步骤是图像缩放，我也将图像缩放实现了cuda并行化操作。下面以这个简单的例子谈谈cuda并行化编程的感受。
         首先肯定是要了解图像缩放的算法，我参考了下面的资料，了解了双线性插值算法。
双线性插值实现的缩放 http://blog.csdn.net/qiqi5521/article/details/2207562 http://www.cnblogs.com/funny-world/p/3162003.html
 
       简单地说，其实就是，定义原图是src,目标图像是dst。遍历目标图像的每个像素的位置（x,y）,按照一定的比例从原图src中找该点（x,y）的附近的四个点，然后根据一定的权重和四个点的像素值算出目标图像（x,y）的像素值。权重的计算就是双线性插值了。
 
      那么怎么并行化这个算法呢？显然，每个位置像素的计算都是独立的，互补影响，就可以一个线程计算一个像素了。不过从实际编程中，我发现有一个地方会严重影响性能。那就是，线程块的数量和每个线程块里面线程的数量。不同的设定，性能很不同。我测试了很多种设定，找到了最优的。但是找最优的设定并没有规律可寻，貌似只能实验。
 
       然后我开始搜索相关资料，看了下面这个 http://stackoverflow.com/questions/11592450/how-to-adjust-the-cuda-number-of-block-and-of-thread-to-get-optimal-performances Howto adjust the cuda number of block and of thread to
get optimalperformances
 
       并没有找到我想要的答案。不过知道了调整线程块的数量和每个线程块线程的数量是性能调优的一方面。我于是又搜索了一下，发现有不少性能调优的文档。带着编程问题去搜索资料这种学习方式，看来行得通。不过cuda里面的官方文档还是要看，特别是programmingguide和bestpractice。
 
       下面是我自己实现的cuda缩放图像的代码：

#include "../../common/book.h"
#include "../../common/cuPrintf.cu"
//HandleError
#include "cuImage.h"

__global__ void resizeGPU(const unsigned char*src,int srcWidth,int srcHeight,unsigned char *dst,int dstWidth,int dstHeight)
{
double srcXf;
double srcYf;
int srcX;
int srcY;
double u;
double v;
int dstOffset;

int y = blockIdx.y*blockDim.y+threadIdx.y;
int x = blockIdx.x*blockDim.x+threadIdx.x;

if(x>=dstWidth || y>=dstHeight) return;

srcXf=  x* ((float)srcWidth/dstWidth) ;//这个float重要阿，不然会算是精度的，尼玛
srcYf =  y* ((float)srcHeight/dstHeight);
srcX = (int)srcXf;
srcY = (int)srcYf;
u= srcXf - srcX;
v = srcYf - srcY;

//r chanel
dstOffset =(y*dstWidth+x)*3;
dst[dstOffset] = 0;
dst[dstOffset]+=(1-u)*(1-v)*src[(srcY*srcWidth+srcX)*3];
dst[dstOffset]+=(1-u)*v*src[((srcY+1)*srcWidth+srcX)*3];
dst[dstOffset]+=u*(1-v)*src[(srcY*srcWidth+srcX+1)*3];
dst[dstOffset]+= u*v*src[((srcY+1)*srcWidth+srcX+1)*3];

//g chanel
dstOffset =(y*dstWidth+x)*3+1;
dst[dstOffset] = 0;
dst[dstOffset]+=(1-u)*(1-v)*src[(srcY*srcWidth+srcX)*3+1];
dst[dstOffset]+=(1-u)*v*src[((srcY+1)*srcWidth+srcX)*3+1];
dst[dstOffset]+=u*(1-v)*src[(srcY*srcWidth+srcX+1)*3+1];
dst[dstOffset]+= u*v*src[((srcY+1)*srcWidth+srcX+1)*3+1];

//b chanel
dstOffset =(y*dstWidth+x)*3+2;
dst[dstOffset] = 0;
dst[dstOffset]+=(1-u)*(1-v)*src[(srcY*srcWidth+srcX)*3+2];
dst[dstOffset]+=(1-u)*v*src[((srcY+1)*srcWidth+srcX)*3+2];
dst[dstOffset]+=u*(1-v)*src[(srcY*srcWidth+srcX+1)*3+2];
dst[dstOffset]+= u*v*src[((srcY+1)*srcWidth+srcX+1)*3+2];

}
void cuResize(const unsigned char*src,int srcWidth,int srcHeight,unsigned char *dst,int dstWidth,int dstHeight)
{

// capture the start time
cudaEvent_t     start, stop;
HANDLE_ERROR( cudaEventCreate( &start ) );
HANDLE_ERROR( cudaEventCreate( &stop ) );
HANDLE_ERROR( cudaEventRecord( start, 0 ) );

int uint = 16;//测了很多次，发现这个比例是最佳的
dim3 grid((dstWidth+uint-1)/uint,(dstHeight+uint-1)/uint);
dim3 block(uint,uint);

unsigned char *dev_src;
HANDLE_ERROR(cudaMalloc((void**)&dev_src,srcWidth*srcHeight*3*sizeof(unsigned char)));
HANDLE_ERROR(cudaMemcpy(dev_src,src,srcWidth*srcHeight*3*sizeof(unsigned char),cudaMemcpyHostToDevice));
unsigned char *dev_dst;
HANDLE_ERROR(cudaMalloc((void**)&dev_dst,dstWidth*dstHeight*3*sizeof(unsigned char)));

resizeGPU<<<grid,block>>>(dev_src, srcWidth, srcHeight,dev_dst, dstWidth, dstHeight);
HANDLE_ERROR(cudaMemcpy(dst,dev_dst,dstWidth*dstHeight*3*sizeof(unsigned char),cudaMemcpyDeviceToHost));

// get stop time, and display the timing results
HANDLE_ERROR( cudaEventRecord( stop, 0 ) );
HANDLE_ERROR( cudaEventSynchronize( stop ) );
float   elapsedTime;
HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime,
start, stop ) );
printf( "GPU:Time to generate:  %3.1f ms\n", elapsedTime );

HANDLE_ERROR( cudaEventDestroy( start ) );
HANDLE_ERROR( cudaEventDestroy( stop ) );

HANDLE_ERROR(cudaFree( dev_src ));
HANDLE_ERROR(cudaFree( dev_dst ));

}