cuda版本的meanshift算法，其核函数执行效率太低

cuda版本的meanshift算法，其核函数执行效率太低

« 于: 七月 03, 2018, 04:09:54 pm »

大家好！我正在实现一个cuda版本的能处理4通道数据的mean shift算法。算法处理1920x1080的4通道图像。跑下来居然需要几十秒钟。opencv的meanShiftProc实现相似的功能，只需要1秒左右吧。大家帮我看看我的代码，帮我分析分析，谢谢！

程序代码: [选择]

__global__ void meanshift_filter_kernel(char *range, char *spatial, int width, int height, size_t rpitch, size_t spitch, int sr, int sp, int maxiter, float eps) { int x0; int y0; x0 = blockIdx.x * blockDim.x + threadIdx.x; y0 = blockIdx.y * blockDim.y + threadIdx.y; if (x0 < width && y0 < height) { do_meanshift(x0, y0, range, spatial, rpitch, spitch, width, height, sr, sp, maxiter, eps); } } __device__ void do_meanshift(int x0, int y0, char *range, char *spatial, int width, int height, size_t rpitch, size_t spitch, int sr, int sp, int maxiter, float eps) { int sr2; int top; int bottom; int left; int right; uchar4 r0; uchar4 r1; int dist; int rmx; int rmy; int rmz; int rmw; int spmx; int spmy; int total; float k; int flag; int x; int y; int i; int x00; int y00; sr2 = sr * sr; flag = 0; r0 = tex2D(intext, x0, y0); x00 = x0; y00 = y0; for (i = 0; i < maxiter; i++) { top = y0 - sp; bottom = y0 + sp; left = x0 - sp; right = x0 + sp; rmx = 0; rmy = 0; rmz = 0; rmw = 0; spmx = 0; spmy = 0; total = 0; for (y = top; y <= bottom; y++) { for (x = left; x <= right; x++) { r1 = tex2D(intext, x, y); dist = (r0.x - r1.x) * (r0.x - r1.x) + (r0.y - r1.y) * (r0.y - r1.y) + \ (r0.z - r1.z) * (r0.z - r1.z) + (r0.w - r1.w) * (r0.w - r1.w); if (dist <= sr2) { rmx += r1.x; rmy += r1.y; rmz += r1.z; rmw += r1.w; spmx += x; spmy += y; total++; } } } if (0 == total) { break; } k = 1.f / total; rmx = __float2int_rz(rmx * k); rmy = __float2int_rz(rmy * k); rmz = __float2int_rz(rmz * k); rmw = __float2int_rz(rmw * k); spmx = __float2int_rz(spmx * k); spmy = __float2int_rz(spmy * k); dist = (rmx - r0.x) * (rmx - r0.x) + (rmy - r0.y) * (rmy - r0.y) + \ (rmz - r0.z) * (rmz - r0.z) + (rmw - r0.w) * (rmw - r0.w); flag = ((spmx == x0 && spmy == y0) || \ (abs(spmx - x0) + abs(spmy - y0) + dist <= eps)); r0.x = rmx; r0.y = rmy; r0.z = rmz; r0.w = rmw; x0 = spmx; y0 = spmy; if (flag) { break; } } *((uchar4 *)(range + y00 * rpitch * sizeof(uchar4) + x00 * sizeof(uchar4))) = r0; *((short2 *)(spatial + y00 * spitch + x00 * sizeof(short2))) = make_short2(x0, y0); }

已记录

« 回复 #1 于: 七月 04, 2018, 01:01:30 pm »

Hi, 几道,

你的代码经过一天时间, 已经阅读完毕.

有几个方面可能存在问题:

(1)根据你的文字描述: "算法处理1920x1080的4通道图像". 但是根据你的profiler的截图中泄漏的信息(你并没有直接告诉我们你的kernel的启动形状配置). 你只使用了(32,8)的block形状, 和(8, 32)的grid形状, 而你的中心点计算方式为:
x0 = blockIdx.x * blockDim.x + threadIdx.x
和: y0 = blockIdx.y * blockDim.y + threadIdx.y;
这将只能覆盖256x256的区域, 和你的处理1920x1080的4通道图像的说法不符合.

(2)你的代码完全没有注释, 你也没有给出任何说法. 你只给出了算法名字(meanshift), 请至少简单说明一下情况.
此外, 将近两个数量级上的性能差别(你的几十秒 vs OpenCV的1秒), 往往代码你在算法层面可能出现问题. 其次才是实现上. 你应当找到你目前所能找到的算法或者说数学上的最佳方案, 然后才应当考虑实现问题.

(3)当你确定下来最佳算法, 并且写出了实现(假设本文是), 你应当至少看下profiler报告的是你在你的当前的卡上, kernel运行起来是访存密集型, 还是计算密集. 然后才能有针对的处理. 这也是正常开发的常规流程.

----------------关于你的代码的具体部分-------------------

(4)Maxwell/Pascal(你是前者, 信息来自你的profiler图片)系列卡的浮点数处理是全速率的, 但是32-bit整数不是. 你大量在使用int计算, 性能上是吃亏的. 例如你的代码写在最前面的海量int变量定义.

(5)代码不应当使用老式C89风格, ;你的变量的时候直接给出定义, 而不是在遥远的前方, 例如看到dist = ....;这行的时候, 需要反复往前扒拉看看dist是什么鬼. 结果发现是一个int变量. 这无论对你(代码书写者), 我(阅读代码的人), 你的同事(后来维护代码的人)都是一个灾难.

(6)看到这代码最终生成了使用非常多寄存器的编译结果(你最后的图片中有), 而代码中很多是对4个通道的这样计算:
dist = (r0.x - r1.x) * (r0.x - r1.x) + (r0.y - r1.y) * (r0.y - r1.y) + (r0.z - r1.z) * (r0.z - r1.z) + (r0.w - r1.w) * (r0.w - r1.w);
应当考虑使用显卡自带的视频处理指令(也可以用于处理图像), 例如说:
4个8-bit的通道被存储于1个uint32_t的话(和你的uchar4一样的, 就地免费类型转换一下即可), 分别进行独立求差, 然后求距离(的平方).
可以分别考虑:
(a)4个分量的差的绝对值: uint32_t c = __vabsdiffu4(a,b);
将计算a和b的4个通道单独相减, 然后进行绝对值操作. 这样比你4次分量提取, 转换成有符号int, 减法, 要好的多. 节省了大量的寄存器和指令.
(b)点积(用于距离的平方)
uint32_t distance = __dp4a(c,c, 0); 将c的4个分量(之前的a和b的4个分量的差的绝对值)进行平方, 然后累加.

这样的写法远比你的原始写法节省寄存器和效率高.

(7) 不应当使用过老的卡, 例如你的GTX750 Ti(信息来自你的Profiler截图). 这是一张目前只有100-200元的卡, 一般做我们这行的, 如果上不起几万元级别的卡(例如Titan-V), 至少也应当买个几千元的卡, 这是常见情况. 只使用100多元(你没看错0的个数)的卡, 连回复本主题的人力成本都不够. 甚至刚才(6)点中的建议的第二点, 就需要普通的6.1的Pascal, 例如GTX1080/1070/1060/1050/1030, 你的750是不能支持的(你的卡可以依然使用前一部分的绝对差). 这也是对你的时间的浪费.

(8)CUDA的并行计算, 尽量应当让线程或者warps具有同样的工作强度, 不能其他线程工作少, 都提前结束了, 突然有部分线程具有额外的工作, 在那里突击. 这样的不均衡不会有好性能的. 例如你看你的这里:

程序代码: [选择]

for (int i = 0; i < max_iter; i++) { .....计算出来了total.... if (total == 0) break; .... 和 if (flat) break; }

这个实际上在线程间的进展是不可预估的(依赖于每个线程得到的total和flag的值的情况). 无法确保线程间尽量一致的工作负载. 请尽量不要这样做.

(9)根据你的代码, 你似乎是在一个小的滑动窗口附近搜索什么东西(或者什么条件). 根据你的代码里面的第二重循环, 将在(x0, y0)点附近的一个大约2SP x 2SP的区域(大致, 这里不是很精确)对所有的点进行某种分量空间上的距离(的平方)计算, 然后统计出来有多少个满足条件(小于SR2)的点, 将这些点的所有分量值进行平均 (你的代码的rmx, rmy, rmz, rmw累加后平均), 将坐标也进行平均, 得到一个虚拟新点(坐标来自平均值), 和虚拟分量(同样来自平均值), 然后重新搜索, 直到达到最大iter的值的时候, 或者达到满足新点固定不动了, 或者坐标和dist足够接近了中的任何一个条件就终止搜索.

这个处理过程应当考虑数学(算法)上的优化的. 此外, 使用warp甚至block级别来集体完成这个, 比一个线程独立的串行的重复多次, 也许性能上要好一些. 请思考.

(10) 之前看到有人做过类似算法, 也在本论坛. 大致是使用了高斯金字塔, 进行比例缩放. 因为常规图片临近的颜色是相近的. 可以在多层分别缩减了, 例如64X, 16X, 4X, 的图片上, 分别求出它们之上的最接近的位置, 然后在进一步的细节图片上更加逼近一点. 也许会显著的减少运算量. 而不是你这样上去在较大的原图上直接搜索? 注意我不是图像处理专业的, 这个仅供参考. 你应当阅读相关算法(meanshift)的文章, 选择最佳算法, 然后进行实现.

(11)根据你的代码, 窗口的大小(SP)变化会显著的加大或者减少运算量. 请确保你的代码和用来参考的实现使用了相同的参数(OpenCV的1s).

目前给出11条建议. 因为你的论坛标志表明你并未购买我们的任何产品或者服务. 因此只能给出这些.

Regards,
屠戮人神.

已记录

« 回复 #2 于: 七月 05, 2018, 12:35:45 pm »

谢谢[名词2]的宝贵建议！对于[名词2]提到：（1）均值偏移（mean shift filter）算法的目标是能实时处理1920x1080分辨率的4通道图像（RGBA或者其他特征图），目前代码速度过慢，所以贴出来的结果是处理256x256图像的结果。

现在我作了如下修改：
（1）值域图（range）由unsigned char型更改为uchar4型，收敛点/空域图（spatial）由unsigned char型改为short2型（存储x和y坐标）。
（2）取消迭代中断，所有线程统一迭代maxiter次。其他两个条件语句保留。
（3）其他修改。

这样跑下来，divergent branch警告得以消除，warp execution efficiency得以提升，但global memory store仍然很低。难道global memory的访问没有实现
coalesced？

以下是修改后的代码：
------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

程序代码: [选择]

__global__ void meanshift_filter_kernel(uchar4 *range, short2 *spatial, int width, int height, int sr, int sp, int maxiter, float eps) { int x0 = blockIdx.x * blockDim.x + threadIdx.x; int y0 = blockIdx.y * blockDim.y + threadIdx.y; if (x0 < width && y0 < height) { do_meanshift(x0, y0, range, spatial, width, height, sr, sp, maxiter, eps); } } __device__ void do_meanshift(int x0, int y0, uchar4 *range, short2 *spatial, int width, int height, int sr, int sp, int maxiter, float eps) { int sr2 = sr * sr; uchar4 r0 = tex2D(intext, x0, y0); for (int i = 0; i < maxiter; i++) { int top = y0 - sp; int bottom = y0 + sp; int left = x0 - sp; int right = x0 + sp; int rmx = 0; int rmy = 0; int rmz = 0; // int rmw = 0; int spmx = 0; int spmy = 0; int total = 0; for (int y = top; y <= bottom; y++) { for (int x = left; x <= right; x++) { uchar4 r1 = tex2D(intext, x, y); int dist = (r0.x - r1.x) * (r0.x - r1.x) + (r0.y - r1.y) * (r0.y - r1.y) + \ (r0.z - r1.z) * (r0.z - r1.z); // +(r0.w - r1.w) * (r0.w - r1.w); if (dist <= sr2) { rmx += r1.x; rmy += r1.y; rmz += r1.z; //rmw += r1.w; spmx += x; spmy += y; total++; } //bool flag = (r0.x - r1.x) * (r0.x - r1.x) + (r0.y - r1.y) * (r0.y - r1.y) + \ // (r0.z - r1.z) * (r0.z - r1.z)/* + (r0.w - r1.w) * (r0.w - r1.w)*/ <= sr2; // //rmx += r1.x * flag; //rmy += r1.y * flag; //rmz += r1.z * flag; rmw += r1.w * flag; //spmx += x * flag; //spmy += y * flag; // //total += flag; } } if (0 == total) { break; } //bool flag = 0 == total; //total += flag; float k = 1.f / total; rmx = __float2int_rz(rmx * k); rmy = __float2int_rz(rmy * k); rmz = __float2int_rz(rmz * k); // rmw = __float2int_rz(rmw * k); spmx = __float2int_rz(spmx * k); spmy = __float2int_rz(spmy * k); /*int dist = (rmx - r0.x) * (rmx - r0.x) + (rmy - r0.y) * (rmy - r0.y) + \ (rmz - r0.z) * (rmz - r0.z) + (rmw - r0.w) * (rmw - r0.w); bool stop = ((spmx == x0 && spmy == y0) || \ (::abs(spmx - x0) + ::abs(spmy - y0) + dist <= eps));*/ r0.x = rmx; r0.y = rmy; r0.z = rmz; //r0.w = rmw; x0 = spmx; y0 = spmy; //int reverse = !flag; // //r0.x = reverse * rmx + flag * r0.x; //r0.y = reverse * rmy + flag * r0.y; //r0.z = reverse * rmz + flag * r0.z; r0.w = reverse * rmw + flag * r0.w; // //x0 = reverse * spmx + flag * x0; //y0 = reverse * spmy + flag * y0; //if (stop) { // break; //} } *(range + (blockIdx.y * blockDim.y + threadIdx.y) * width + \ (blockIdx.x * blockDim.x + threadIdx.x)) = r0; *(spatial + (blockIdx.y * blockDim.y + threadIdx.y) * width + \ (blockIdx.x * blockDim.x + threadIdx.x)) = make_short2(x0, y0); } 作为对比，我把第四通道（w分量）的处理屏蔽了，并贴出OpenCV的实现： __global__ void meanshiftproc_kernel(unsigned char* outr, size_t outrstep, unsigned char* outsp, size_t outspstep, int cols, int rows, int sp, int sr, int maxIter, float eps) { int x0 = blockIdx.x * blockDim.x + threadIdx.x; int y0 = blockIdx.y * blockDim.y + threadIdx.y; if( x0 < cols && y0 < rows ) { int basesp = (blockIdx.y * blockDim.y + threadIdx.y) * outspstep + (blockIdx.x * blockDim.x + threadIdx.x) * 2 * sizeof(short); *(short2*)(outsp + basesp) = do_mean_shift(x0, y0, outr, outrstep, cols, rows, sp, sr, maxIter, eps); } } __device__ short2 do_mean_shift(int x0, int y0, unsigned char* out, size_t out_step, int cols, int rows, int sp, int sr, int maxIter, float eps) { int isr2 = sr*sr; uchar4 c = tex2D(tex_meanshift, x0, y0 ); // iterate meanshift procedure for( int iter = 0; iter < maxIter; iter++ ) { int count = 0; int s0 = 0, s1 = 0, s2 = 0, sx = 0, sy = 0; float icount; //mean shift: process pixels in window (p-sigmaSp)x(p+sigmaSp) int minx = x0-sp; int miny = y0-sp; int maxx = x0+sp; int maxy = y0+sp; for( int y = miny; y <= maxy; y++) { int rowCount = 0; for( int x = minx; x <= maxx; x++ ) { uchar4 t = tex2D( tex_meanshift, x, y ); int norm2 = (t.x - c.x) * (t.x - c.x) + (t.y - c.y) * (t.y - c.y) + (t.z - c.z) * (t.z - c.z); if( norm2 <= isr2 ) { s0 += t.x; s1 += t.y; s2 += t.z; sx += x; rowCount++; } } count += rowCount; sy += y*rowCount; } if( count == 0 ) break; icount = 1.f/count; int x1 = __float2int_rz(sx*icount); int y1 = __float2int_rz(sy*icount); s0 = __float2int_rz(s0*icount); s1 = __float2int_rz(s1*icount); s2 = __float2int_rz(s2*icount); int norm2 = (s0 - c.x) * (s0 - c.x) + (s1 - c.y) * (s1 - c.y) + (s2 - c.z) * (s2 - c.z); bool stopFlag = (x0 == x1 && y0 == y1) || (::abs(x1-x0) + ::abs(y1-y0) + norm2 <= eps); x0 = x1; y0 = y1; c.x = s0; c.y = s1; c.z = s2; if( stopFlag ) break; } int base = (blockIdx.y * blockDim.y + threadIdx.y) * out_step + (blockIdx.x * blockDim.x + threadIdx.x) * 4 * sizeof(uchar); *(uchar4*)(out + base) = c; return make_short2((short)x0, (short)y0); }

您可以看到，这两份代码非常接近，然而运行效率却相差太大。下面是profiler分析结果。

我的代码：

OpenCV的代码：

已记录

« 回复 #3 于: 七月 05, 2018, 03:00:21 pm »

不想看太多了，反正我写的你也不看。

如果你认为OpenCV的开源代码就是最优实现，你完全可以直接无视我之前的任何建议，而采用OpenCV的实现。反正“代码是一样的“，人家却能运行快100倍。果断直接抄。

就这样了。

已记录

« 回复 #4 于: 七月 05, 2018, 03:04:14 pm »

但需要额外说明的是，这种依赖于输入的算法，运行最终时间可能不定（和你的输入有关），

例如在快速能达到count == 0的时候，和你循环了很多次才能达到的时候。

不妨确保在相同的参数下，包括图像，范围参数，编译参数等等，这些都实验后再考虑。但是这些原因应当是你自己去负责的，而不是我在这里猜测。

我不能继续看你的代码。抱歉。

已记录

« 回复 #5 于: 七月 05, 2018, 04:58:53 pm »

本帖最后由几道于 2018-7-5 17:01 编辑

[名词2]的回复我是认真看了的，也查阅的很多相关的资料:)

OpenCV的实现当然不是最优的，它没有任何搜索策略的优化。定位跑得慢的原因是关键，这样才能逐步提高代码的性能，也才能在后期考虑更多数学上的优化。当我把OpenCV源码之中mean shift部分的代码单独提取出来编译的时候（不调用库中对应的接口），有趣的是，其运行速度比直接调用相应的接口慢100来倍，和我自己的代码差不多！这个有点让人疑惑了。我把OpenCV编译mean shift模块时的NVCC的配置文件（cmake）打开看了看，也没看到什么特别的差异——和我自己的工程配置。

不管怎么样，谢谢[名词2]的耐心解答。

已记录

引用自: 几道于七月 05, 2018, 04:58:53 pm

« 回复 #6 于: 七月 05, 2018, 06:23:03 pm »

[名词2]的回复我是认真看了的，也查阅的很多相关的资料

OpenCV的实现当然不是最优的，它没有任何搜索策 ...

再重复一次，

你确定输入参数（图像），方形范围大小，编译参数（例如优化，计算能力设定），显卡，这些都一样么？

这个实际上上面已经问过一次了，但是依然被你照常的无视了。。。。

已记录