斑竹你好,在写图像相关的代码中遇到个问题。代码如下:
__global__ void CorrelationGPU(float *Img1 ,float *Img2,int Height,int Width,int *SubMatrix,int M,int m_nRow,int m_nCol)
{
int bky=blockIdx.y;
int bkx=blockIdx.x;
int tdy=threadIdx.y;
int tdx=threadIdx.x;
int i,j;
int y=SubMatrix[bky*m_nRow+bkx+m_nRow*m_nCol]; //原图block的中心点
int x=SubMatrix[bky*m_nRow+bkx];
__shared__ float x1,x2,x3;
__shared__ float shared6[4847]; //大小是(M/2*2+1)*(M/2*2+1)*3,存储两个模板的值
M=M/2;
i=tdy-M;
j=tdx-M; //以块的中心为坐标
int M2=(2*M+1)*(2*M+1); //让块取成奇数*奇数,方便用
//////////////////////////////////////////////////////////////////////////
//求以(x,y)为中心的,(2M+1)*(2M+1)为模板的的相关系数
shared6[0*M2+tdy*blockDim.x+tdx]=Img1[(y+i)*Width+x+j]; //f
shared6[1*M2+tdy*blockDim.x+tdx]=Img2[(y+i)*Width+x+j]; //g
shared6[2*M2+tdy*blockDim.x+tdx]=shared6[0*M2+tdy*blockDim.x+tdx]; //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=shared6[1*M2+tdy*blockDim.x+tdx];
__syncthreads();
for (i=M2/2;i>0;i=i>>1){ //奇数做规约,这里用i节省空间
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0];
x2=shared6[3*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1]; //做的奇数规约x1~fm,x2~gm
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x1=x1/(M2);
x2=x2/(M2);
//用规约求相关系数cor
shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[0*M2+tdy*blockDim.x+tdx]-x1); //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-x2)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
__syncthreads();
for (i=M2/2;i>0;i=i>>1){//规约,得到fsum,gsum,fgsum
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
shared6[4*M2+tdy*blockDim.x+tdx]+=shared6[4*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0]; //x1~fsum/cor,x2~gsum,x3~fgsum
x2=shared6[3*M2+0];
x3=shared6[4*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1];
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x3+=shared6[4*M2+blockDim.y*blockDim.x-1];
x1=x3/sqrt(x1*x2);
}
代码说明如下:
1:其中Img1,Img2是两幅图,高和宽分别为Height、Width,把两个图分割成m_nRowm_nCol份,每份的大小是MM,中心点存储在数组SubMatrix中,这里M=31。
2:在核函数中申请了一个共享数组,用来存储Img1、Img2的数据,以便用来计算相关系数;
3:因为线程是块的左上角为原点的,为了方便这里用i=ty-M;j=tx-M,把坐标原点放到了块的中心位置;
4:做相关的公式如下:
[attach]3248[/attach]
具体描述在原来这个帖子中有:http://cudazone.nvidia.cn/forum/forum.php?mod=viewthread&tid=7001
5:因此先用规约求均值:
shared6[0*M2+tdy*blockDim.x+tdx]=Img1[(y+i)*Width+x+j]; //f
shared6[1*M2+tdy*blockDim.x+tdx]=Img2[(y+i)*Width+x+j]; //g
shared6[2*M2+tdy*blockDim.x+tdx]=shared6[0*M2+tdy*blockDim.x+tdx]; //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=shared6[1*M2+tdy*blockDim.x+tdx];
__syncthreads();
这里第一部分是传入数据,第二部分是把数据拷贝到共享数据的另一块,因为在后面求相关系数的时候还要用,但是在递归求和的时候要改变数组里面的数据,所以这样是为了保护原有数据;
6:这部分就是用规约求均值:
for (i=M2/2;i>0;i=i>>1){ //奇数做规约,这里用i节省空间
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){ //这里是把和值给所有的线程,如果给share,后面再除的时候要出问题(出现多次除法)
x1=shared6[2*M2+0];
x2=shared6[3*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1]; //做的奇数规约x1~fm,x2~gm
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x1=x1/(M2);
x2=x2/(M2);
因为数组为奇数,所以在后面补加了一项。在这里得到的x1和x2就是两个块的均值。
7:但是在程序这里出了问题:
shared6[2*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[0*M2+tdy*blockDim.x+tdx]-x1); //数据,这里传给0,1,2是为了不让规约运算改变原始值
shared6[3*M2+tdy*blockDim.x+tdx]=(shared6[1*M2+tdy*blockDim.x+tdx]-x2)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
shared6[4*M2+tdy*blockDim.x+tdx]=(shared6[0*M2+tdy*blockDim.x+tdx]-x1)*(shared6[1*M2+tdy*blockDim.x+tdx]-x2);
__syncthreads();
for (i=M2/2;i>0;i=i>>1){//规约,得到fsum,gsum,fgsum
if (tdy*blockDim.x+tdx<i){
shared6[2*M2+tdy*blockDim.x+tdx]+=shared6[2*M2+tdy*blockDim.x+tdx+i];
shared6[3*M2+tdy*blockDim.x+tdx]+=shared6[3*M2+tdy*blockDim.x+tdx+i];
shared6[4*M2+tdy*blockDim.x+tdx]+=shared6[4*M2+tdy*blockDim.x+tdx+i];
}
__syncthreads();
}
if (tdy*blockDim.x+tdx==0){
x1=shared6[2*M2+0]; //x1~fsum/cor,x2~gsum,x3~fgsum
x2=shared6[3*M2+0];
x3=shared6[4*M2+0];
}
x1+=shared6[2*M2+blockDim.y*blockDim.x-1];
x2+=shared6[3*M2+blockDim.y*blockDim.x-1];
x3+=shared6[4*M2+blockDim.y*blockDim.x-1];
x1=x3/sqrt(x1*x2);
在第6步中x1,x2的值还是正确的,这里是58.064,但是在往下面走,到第7步中的__syncthreads值就变成了0.11869,这里明显出现错误。但是如果把__shared__ float x1,x2,x3定义成float x1,x2,x3就不会出现上面的错误。
8:我想是不是在这里因为x1,x2,x3是共享的,所以x1=x1/(M2);x2=x2/(M2) 被执行了多次,但是只是我猜测而已,还请版主大人指教。
9:我这里要用共享空间的原因是往后面做我是要用上面的数据得到一个方程的解,比如AX=B,其中A、B都是共享空间上的,A是66的,B是61的,通过求解出来的X值判断一个块内循环是否继续。即所有块内的线程就以这个共享的空间x为判断标准。(这个可能没见多明白,但是感觉和问题关系不大。。谢谢)