问题描述
|
这是我在CUDA中实现的伪代码的一部分,作为图像重建算法的一部分:
for each xbin(0->detectorXDim/2-1):
for each ybin(0->detectorYDim-1):
rayInit=(xbin*xBinSize+0.5,ybin*xBinSize+0.5,-detectordistance)
rayEnd=beamFocusCoord
slopeVector=rayEnd-rayInit
//kNowing that r=rayInit+t*slopeVector;
//x=rayInit[0]+t*slopeVector[0]
//y=rayInit[1]+t*slopeVector[1]
//z=rayInit[2]+t*slopeVector[2]
//to find ray xx intersections:
for each xinteger(xbin+1->detectorXDim/2):
solve t for x=xinteger*xBinSize;
find corresponding y and z
add to intersections array
//find ray yy intersections(analogous to xx intersections)
//find ray zz intersections(analogous to xx intersections)
到目前为止,这是我想出的:
__global__ void sysmat(int xfocus,int yfocus,int zfocus,int xbin,int xbinsize,int ybin,int ybinsize,int zbin,int projecoes){
int tx=threadIdx.x,ty=threadIdx.y,tz=threadIdx.z,bx=blockIdx.x,by=blockIdx.y,i,x,y,z;
int idx=ty+by*blocksize;
int idy=tx+bx*blocksize;
int slopeVectorx=xfocus-idx*xbinsize+0.5;
int slopeVectory=yfocus-idy*ybinsize+0.5;
int slopeVectorz=zfocus-zdetector;
__syncthreads();
//points where the ray intersects x axis
int xint=idx+1;
int yint=idy+1;
int*intersectionsx[(detectorXDim/2-xint)+(detectorYDim-yint)+(zfocus)];
int*intersectionsy[(detectorXDim/2-xint)+(detectorYDim-yint)+(zfocus)];
int*intersectionsz[(detectorXDim/2-xint)+(detectorYDim-yint)+(zfocus)];
for(xint=xint; xint<detectorXDim/2;xint++){
x=xint*xbinsize;
t=(x-idx)/slopeVectorx;
y=idy+t*slopeVectory;
z=z+t*slopeVectorz;
intersectionsx[xint-1]=x;
intersectionsy[xint-1]=y;
intersectionsz[xint-1]=z;
__syncthreads();
}
...
}
这只是一部分代码。我知道可能会有一些错误(如果它们很明显是错误的,您可以指出这些错误),但是我更担心的是:
每个线程(对应于检测器仓)需要三个阵列,因此它可以保存射线(穿过此线程/仓)与x,y和z轴的倍数相交的点。每个阵列的长度取决于检测器中线程/容器(其索引)的位置以及beamFocusCoord(固定)。为了做到这一点,我编写了这段代码,我确定这是不可能完成的(用一个小的测试内核确认了它,并返回错误:\“表达式必须具有恒定值\”):
int*intersectionsx[(detectorXDim/2-xint)+(detectorXDim-yint)+(zfocus)];
int*intersectionsy[(detectorXDim/2-xint)+(detectorXDim-yint)+(zfocus)];
int*intersectionsz[(detectorXDim/2-xint)+(detectorXDim-yint)+(zfocus)];
所以最后,我想知道这段代码是否有替代方法,其中向量的长度取决于分配该向量的线程的索引。
先感谢您 ;)
编辑:鉴于每个线程将必须保存一个数组,其中包含射线(从光束源到检测器)与xx,yy和zz轴之间的交点的坐标,并且空间尺寸在(I目前没有确切的数字,但它们非常接近实际值(1400x3600x60),使用CUDA可以解决此问题吗?
例如,线程(0,0)在x轴上具有1400个交集,在y轴上具有3600个交集,在z轴上具有60个交集,这意味着我将不得不创建一个大小为(1400 + 3600 + 60)*的数组sizeof(float),每个线程大约20kb。
因此,鉴于每个线程都超过16kb的本地内存,这是不可能的。另一种选择是分配这些数组,但是通过更多的数学运算,我们得到(1400 + 3600 + 60)* 4 *线程数(即1400 * 3600),这也超过了可用的全局内存量:(
因此,我没有足够的想法来解决这个问题,我们将为您提供任何帮助。
解决方法
没有。
必须在内核启动时知道CUDA中的每个内存。内核运行时,您无法分配/取消分配/更改任何内容。对于全局存储器,共享存储器和寄存器,这是正确的。
常见的解决方法是预先分配所需的最大内存大小。这可以简单地多次分配一个线程所需的最大大小,也可以复杂到将所有那些线程需要的大小相加以获得最大的总和并在该数组中计算适当的线程偏移量。这是内存分配和偏移量计算时间之间的折衷方案。
由于内存有限,请尽量选择简单的解决方案,而必须这样做则选择复杂的解决方案。
,为什么不使用纹理?使用2D或3D纹理将使此问题容易得多。 GPU旨在执行非常快速的浮点插值,并且CUDA对此提供了出色的支持。文献中有在GPU上进行投影重建的示例,例如使用支持CUDA的GPU加速具有运动补偿的同时代数重构技术,纹理是其算法不可或缺的一部分。除非您需要像正弦插值之类的怪异东西,否则您自己的手动坐标计算只会比GPU提供的速度慢且容易出错。
1400x3600x60对于单个3D纹理而言有点大,但是您可以将问题分解为2D切片,3D子体积或分层多分辨率重建。这些都已被其他研究人员使用。只需搜索PubMed。