From 21ad8e92d4b039817e4004bfc952be0576eb4d1d Mon Sep 17 00:00:00 2001 From: vbystricky Date: Mon, 29 Sep 2014 12:46:56 +0400 Subject: [PATCH] Optimize OpenCL version of StereoBM function Fix problems on NVidia devices. --- modules/calib3d/src/opencl/stereobm.cl | 267 ++++++++++++++----------- modules/calib3d/src/stereobm.cpp | 33 ++- 2 files changed, 167 insertions(+), 133 deletions(-) diff --git a/modules/calib3d/src/opencl/stereobm.cl b/modules/calib3d/src/opencl/stereobm.cl index 73402a6a1..86e4c24d3 100644 --- a/modules/calib3d/src/opencl/stereobm.cl +++ b/modules/calib3d/src/opencl/stereobm.cl @@ -44,217 +44,251 @@ ////////////////////////////////////////// stereoBM ////////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// -#ifdef csize - #define MAX_VAL 32767 -void calcDisp(__local short * cost, __global short * disp, int uniquenessRatio, int mindisp, int ndisp, int w, - __local int * bestDisp, __local int * bestCost, int d, int x, int y, int cols, int rows, int wsz2) +#ifndef WSZ +#define WSZ 2 +#endif + +#define WSZ2 (WSZ / 2) + +#ifdef DEFINE_KERNEL_STEREOBM + +#define DISPARITY_SHIFT 4 +#define FILTERED ((MIN_DISP - 1) << DISPARITY_SHIFT) + +void calcDisp(__local short * cost, __global short * disp, int uniquenessRatio, + __local int * bestDisp, __local int * bestCost, int d, int x, int y, int cols, int rows) { - short FILTERED = (mindisp - 1)<<4; - int best_disp = *bestDisp, best_cost = *bestCost, best_disp_back = ndisp - best_disp - 1; + int best_disp = *bestDisp, best_cost = *bestCost; + barrier(CLK_LOCAL_MEM_FENCE); short c = cost[0]; + int thresh = best_cost + (best_cost * uniquenessRatio / 100); + bool notUniq = ( (c <= thresh) && (d < (best_disp - 1) || d > (best_disp + 1) ) ); - int thresh = best_cost + (best_cost * uniquenessRatio/100); - bool notUniq = ( (c <= thresh) && (d < (best_disp_back - 1) || d > (best_disp_back + 1) ) ); - - if(notUniq) + if (notUniq) *bestCost = FILTERED; barrier(CLK_LOCAL_MEM_FENCE); - if( *bestCost != FILTERED && x < cols-wsz2-mindisp && y < rows-wsz2 && d == best_disp_back) + if( *bestCost != FILTERED && x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2 && d == best_disp) { - int y3 = (best_disp_back > 0) ? cost[-w] : cost[w], - y2 = c, - y1 = (best_disp_back < ndisp-1) ? cost[w] : cost[-w]; - int d_aprox = y3+y1-2*y2 + abs(y3-y1); - disp[0] = (short)(((best_disp_back + mindisp)*256 + (d_aprox != 0 ? (y3-y1)*256/d_aprox : 0) + 15) >> 4); + int d_aprox = 0; + int yp =0, yn = 0; + if ((0 < best_disp) && (best_disp < NUM_DISP - 1)) + { + yp = cost[-2 * BLOCK_SIZE_Y]; + yn = cost[2 * BLOCK_SIZE_Y]; + d_aprox = yp + yn - 2 * c + abs(yp - yn); + } + disp[0] = (short)(((best_disp + MIN_DISP)*256 + (d_aprox != 0 ? (yp - yn) * 256 / d_aprox : 0) + 15) >> 4); } } -int calcLocalIdx(int x, int y, int d, int w) -{ - return d*2*w + (w - 1 - y + x); -} - -void calcNewCoordinates(int * x, int * y, int nthread) -{ - int oldX = *x - (1-nthread), oldY = *y; - *x = (oldX == oldY) ? (0*nthread + (oldX + 2)*(1-nthread) ) : (oldX+1)*(1-nthread) + (oldX+1)*nthread; - *y = (oldX == oldY) ? (0*(1-nthread) + (oldY + 1)*nthread) : oldY + 1*(1-nthread); -} - short calcCostBorder(__global const uchar * leftptr, __global const uchar * rightptr, int x, int y, int nthread, - int wsz2, short * costbuf, int * h, int cols, int d, short cost, int winsize) + short * costbuf, int *h, int cols, int d, short cost) { - int head = (*h)%wsz; + int head = (*h) % WSZ; __global const uchar * left, * right; - int idx = mad24(y+wsz2*(2*nthread-1), cols, x+wsz2*(1-2*nthread)); + int idx = mad24(y + WSZ2 * (2 * nthread - 1), cols, x + WSZ2 * (1 - 2 * nthread)); left = leftptr + idx; right = rightptr + (idx - d); - int shift = 1*nthread + cols*(1-nthread); short costdiff = 0; - for(int i = 0; i < winsize; i++) + if (0 == nthread) { - costdiff += abs( left[0] - right[0] ); - left += shift; - right += shift; + #pragma unroll + for (int i = 0; i < WSZ; i++) + { + costdiff += abs( left[0] - right[0] ); + left += cols; + right += cols; + } + } + else // (1 == nthread) + { + #pragma unroll + for (int i = 0; i < WSZ; i++) + { + costdiff += abs(left[i] - right[i]); + } } cost += costdiff - costbuf[head]; costbuf[head] = costdiff; - (*h) = (*h)%wsz + 1; + *h = head + 1; return cost; } short calcCostInside(__global const uchar * leftptr, __global const uchar * rightptr, int x, int y, - int wsz2, int cols, int d, short cost_up_left, short cost_up, short cost_left, - int winsize) + int cols, int d, short cost_up_left, short cost_up, short cost_left) { __global const uchar * left, * right; - int idx = mad24(y-wsz2-1, cols, x-wsz2-1); + int idx = mad24(y - WSZ2 - 1, cols, x - WSZ2 - 1); left = leftptr + idx; right = rightptr + (idx - d); - int idx2 = winsize*cols; + int idx2 = WSZ*cols; uchar corrner1 = abs(left[0] - right[0]), - corrner2 = abs(left[winsize] - right[winsize]), + corrner2 = abs(left[WSZ] - right[WSZ]), corrner3 = abs(left[idx2] - right[idx2]), - corrner4 = abs(left[idx2 + winsize] - right[idx2 + winsize]); + corrner4 = abs(left[idx2 + WSZ] - right[idx2 + WSZ]); return cost_up + cost_left - cost_up_left + corrner1 - corrner2 - corrner3 + corrner4; } -__kernel void stereoBM(__global const uchar * leftptr, __global const uchar * rightptr, __global uchar * dispptr, - int disp_step, int disp_offset, int rows, int cols, int mindisp, int ndisp, - int preFilterCap, int textureTreshold, int uniquenessRatio, int sizeX, int sizeY, int winsize) +__kernel void stereoBM(__global const uchar * leftptr, + __global const uchar * rightptr, + __global uchar * dispptr, int disp_step, int disp_offset, + int rows, int cols, // rows, cols of left and right images, not disp + int textureTreshold, int uniquenessRatio) { - int gx = get_global_id(0)*sizeX; - int gy = get_global_id(1)*sizeY; - int lz = get_local_id(2); + int lz = get_local_id(0); + int gx = get_global_id(1) * BLOCK_SIZE_X; + int gy = get_global_id(2) * BLOCK_SIZE_Y; - int nthread = lz/ndisp; - int d = lz%ndisp; - int wsz2 = wsz/2; + int nthread = lz / NUM_DISP; + int disp_idx = lz % NUM_DISP; __global short * disp; __global const uchar * left, * right; - __local short costFunc[csize]; + __local short costFunc[2 * BLOCK_SIZE_Y * NUM_DISP]; + __local short * cost; __local int best_disp[2]; __local int best_cost[2]; best_cost[nthread] = MAX_VAL; - best_disp[nthread] = MAX_VAL; + best_disp[nthread] = -1; barrier(CLK_LOCAL_MEM_FENCE); - short costbuf[wsz]; + short costbuf[WSZ]; int head = 0; - int shiftX = wsz2 + ndisp + mindisp - 1; - int shiftY = wsz2; + int shiftX = WSZ2 + NUM_DISP + MIN_DISP - 1; + int shiftY = WSZ2; int x = gx + shiftX, y = gy + shiftY, lx = 0, ly = 0; - int costIdx = calcLocalIdx(lx, ly, d, sizeY); + int costIdx = disp_idx * 2 * BLOCK_SIZE_Y + (BLOCK_SIZE_Y - 1); cost = costFunc + costIdx; int tempcost = 0; - if(x < cols-wsz2-mindisp && y < rows-wsz2) + if (x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2) { - int shift = 1*nthread + cols*(1-nthread); - for(int i = 0; i < winsize; i++) + if (0 == nthread) { - int idx = mad24(y-wsz2+i*nthread, cols, x-wsz2+i*(1-nthread)); - left = leftptr + idx; - right = rightptr + (idx - d); - short costdiff = 0; - for(int j = 0; j < winsize; j++) + #pragma unroll + for (int i = 0; i < WSZ; i++) { - costdiff += abs( left[0] - right[0] ); - left += shift; - right += shift; + int idx = mad24(y - WSZ2, cols, x - WSZ2 + i); + left = leftptr + idx; + right = rightptr + (idx - disp_idx); + short costdiff = 0; + for(int j = 0; j < WSZ; j++) + { + costdiff += abs( left[0] - right[0] ); + left += cols; + right += cols; + } + costbuf[i] = costdiff; } - if(nthread==1) + } + else // (1 == nthread) + { + #pragma unroll + for (int i = 0; i < WSZ; i++) { + int idx = mad24(y - WSZ2 + i, cols, x - WSZ2); + left = leftptr + idx; + right = rightptr + (idx - disp_idx); + short costdiff = 0; + for (int j = 0; j < WSZ; j++) + { + costdiff += abs( left[j] - right[j]); + } tempcost += costdiff; + costbuf[i] = costdiff; } - costbuf[head] = costdiff; - head++; } } - if(nthread==1) + if (nthread == 1) { cost[0] = tempcost; - atomic_min(best_cost+nthread, tempcost); + atomic_min(best_cost + 1, tempcost); } barrier(CLK_LOCAL_MEM_FENCE); - if(best_cost[1] == tempcost) - atomic_min(best_disp + 1, ndisp - d - 1); + if (best_cost[1] == tempcost) + atomic_max(best_disp + 1, disp_idx); barrier(CLK_LOCAL_MEM_FENCE); - int dispIdx = mad24(gy, disp_step, disp_offset + gx*(int)sizeof(short)); + int dispIdx = mad24(gy, disp_step, mad24((int)sizeof(short), gx, disp_offset)); disp = (__global short *)(dispptr + dispIdx); - calcDisp(cost, disp, uniquenessRatio, mindisp, ndisp, 2*sizeY, - best_disp + 1, best_cost+1, d, x, y, cols, rows, wsz2); + calcDisp(cost, disp, uniquenessRatio, best_disp + 1, best_cost + 1, disp_idx, x, y, cols, rows); barrier(CLK_LOCAL_MEM_FENCE); lx = 1 - nthread; ly = nthread; - for(int i = 0; i < sizeY*sizeX/2; i++) + for (int i = 0; i < BLOCK_SIZE_Y * BLOCK_SIZE_X / 2; i++) { - x = (lx < sizeX) ? gx + shiftX + lx : cols; - y = (ly < sizeY) ? gy + shiftY + ly : rows; + x = (lx < BLOCK_SIZE_X) ? gx + shiftX + lx : cols; + y = (ly < BLOCK_SIZE_Y) ? gy + shiftY + ly : rows; best_cost[nthread] = MAX_VAL; - best_disp[nthread] = MAX_VAL; + best_disp[nthread] = -1; barrier(CLK_LOCAL_MEM_FENCE); - costIdx = calcLocalIdx(lx, ly, d, sizeY); + costIdx = mad24(2 * BLOCK_SIZE_Y, disp_idx, (BLOCK_SIZE_Y - 1 - ly + lx)); + if (0 > costIdx) + costIdx = BLOCK_SIZE_Y - 1; cost = costFunc + costIdx; - - if(x < cols-wsz2-mindisp && y < rows-wsz2 ) + if (x < cols - WSZ2 - MIN_DISP && y < rows - WSZ2) { - tempcost = ( ly*(1-nthread) + lx*nthread == 0 ) ? - calcCostBorder(leftptr, rightptr, x, y, nthread, wsz2, costbuf, &head, cols, d, - cost[2*nthread-1], winsize) : - calcCostInside(leftptr, rightptr, x, y, wsz2, cols, d, - cost[0], cost[1], cost[-1], winsize); + tempcost = (ly * (1 - nthread) + lx * nthread == 0) ? + calcCostBorder(leftptr, rightptr, x, y, nthread, costbuf, &head, cols, disp_idx, cost[2*nthread-1]) : + calcCostInside(leftptr, rightptr, x, y, cols, disp_idx, cost[0], cost[1], cost[-1]); } cost[0] = tempcost; atomic_min(best_cost + nthread, tempcost); barrier(CLK_LOCAL_MEM_FENCE); - if(best_cost[nthread] == tempcost) - atomic_min(best_disp + nthread, ndisp - d - 1); + if (best_cost[nthread] == tempcost) + atomic_max(best_disp + nthread, disp_idx); barrier(CLK_LOCAL_MEM_FENCE); - int dispIdx = mad24(gy+ly, disp_step, disp_offset + (gx+lx)*(int)sizeof(short)); + dispIdx = mad24(gy + ly, disp_step, mad24((int)sizeof(short), (gx + lx), disp_offset)); disp = (__global short *)(dispptr + dispIdx); - calcDisp(cost, disp, uniquenessRatio, mindisp, ndisp, 2*sizeY, - best_disp + nthread, best_cost + nthread, d, x, y, cols, rows, wsz2); + calcDisp(cost, disp, uniquenessRatio, best_disp + nthread, best_cost + nthread, disp_idx, x, y, cols, rows); + barrier(CLK_LOCAL_MEM_FENCE); - calcNewCoordinates(&lx, &ly, nthread); + if (lx + nthread - 1 == ly) + { + lx = (lx + nthread + 1) * (1 - nthread); + ly = (ly + 1) * nthread; + } + else + { + lx += nthread; + ly = ly - nthread + 1; + } } } - -#endif +#endif //DEFINE_KERNEL_STEREOBM ////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////// Norm Prefiler //////////////////////////////////////////// ////////////////////////////////////////////////////////////////////////////////////////////////// __kernel void prefilter_norm(__global unsigned char *input, __global unsigned char *output, - int rows, int cols, int prefilterCap, int winsize, int scale_g, int scale_s) + int rows, int cols, int prefilterCap, int scale_g, int scale_s) { + // prefilterCap in range 1..63, checked in StereoBMImpl::compute + int x = get_global_id(0); int y = get_global_id(1); - int wsz2 = winsize/2; if(x < cols && y < rows) { @@ -262,13 +296,13 @@ __kernel void prefilter_norm(__global unsigned char *input, __global unsigned ch input[y * cols + max(x-1,0)] * 1 + input[ y * cols + x] * 4 + input[y * cols + min(x+1, cols-1)] * 1 + input[min(y+1, rows-1) * cols + x] * 1; int cov2 = 0; - for(int i = -wsz2; i < wsz2+1; i++) - for(int j = -wsz2; j < wsz2+1; j++) + for(int i = -WSZ2; i < WSZ2+1; i++) + for(int j = -WSZ2; j < WSZ2+1; j++) cov2 += input[clamp(y+i, 0, rows-1) * cols + clamp(x+j, 0, cols-1)]; int res = (cov1*scale_g - cov2*scale_s)>>10; - res = min(clamp(res, -prefilterCap, prefilterCap) + prefilterCap, 255); - output[y * cols + x] = res & 0xFF; + res = clamp(res, -prefilterCap, prefilterCap) + prefilterCap; + output[y * cols + x] = res; } } @@ -280,20 +314,21 @@ __kernel void prefilter_norm(__global unsigned char *input, __global unsigned ch __kernel void prefilter_xsobel(__global unsigned char *input, __global unsigned char *output, int rows, int cols, int prefilterCap) { + // prefilterCap in range 1..63, checked in StereoBMImpl::compute int x = get_global_id(0); int y = get_global_id(1); if(x < cols && y < rows) { - output[y * cols + x] = min(prefilterCap, 255) & 0xFF; - } + if (0 < x && !((y == rows-1) & (rows%2==1) ) ) + { + int cov = input[ ((y > 0) ? y-1 : y+1) * cols + (x-1)] * (-1) + input[ ((y > 0) ? y-1 : y+1) * cols + ((x 0 && !((y == rows-1)&(rows%2==1) ) ) - { - int cov = input[ ((y > 0) ? y-1 : y+1) * cols + (x-1)] * (-1) + input[ ((y > 0) ? y-1 : y+1) * cols + ((xSADWindowSize; int wsz2 = wsz/2; - int sizeX = std::max(11, 27 - ocl::Device::getDefault().maxComputeUnits() ), sizeY = sizeX-1, N = ndisp*2; + ocl::Device devDef = ocl::Device::getDefault(); + int sizeX = devDef.isIntel() ? 32 : std::max(11, 27 - devDef.maxComputeUnits()), + sizeY = sizeX - 1, + N = ndisp * 2; - ocl::Kernel k("stereoBM", ocl::calib3d::stereobm_oclsrc, cv::format("-D csize=%d -D wsz=%d", (2*sizeY)*ndisp, wsz) ); + cv::String opt = cv::format("-D DEFINE_KERNEL_STEREOBM -D MIN_DISP=%d -D NUM_DISP=%d" + " -D BLOCK_SIZE_X=%d -D BLOCK_SIZE_Y=%d -D WSZ=%d", + mindisp, ndisp, + sizeX, sizeY, wsz); + ocl::Kernel k("stereoBM", ocl::calib3d::stereobm_oclsrc, opt); if(k.empty()) return false; @@ -753,15 +760,14 @@ static bool ocl_stereobm( InputArray _left, InputArray _right, int cols = left.cols, rows = left.rows; _disp.create(_left.size(), CV_16S); - _disp.setTo((mindisp - 1)<<4); + _disp.setTo((mindisp - 1) << 4); Rect roi = Rect(Point(wsz2 + mindisp + ndisp - 1, wsz2), Point(cols-wsz2-mindisp, rows-wsz2) ); UMat disp = (_disp.getUMat())(roi); - int globalX = disp.cols/sizeX, globalY = disp.rows/sizeY; - globalX += (disp.cols%sizeX) > 0 ? 1 : 0; - globalY += (disp.rows%sizeY) > 0 ? 1 : 0; - size_t globalThreads[3] = { globalX, globalY, N}; - size_t localThreads[3] = {1, 1, N}; + int globalX = (disp.cols + sizeX - 1) / sizeX, + globalY = (disp.rows + sizeY - 1) / sizeY; + size_t globalThreads[3] = {N, globalX, globalY}; + size_t localThreads[3] = {N, 1, 1}; int idx = 0; idx = k.set(idx, ocl::KernelArg::PtrReadOnly(left)); @@ -769,15 +775,8 @@ static bool ocl_stereobm( InputArray _left, InputArray _right, idx = k.set(idx, ocl::KernelArg::WriteOnlyNoSize(disp)); idx = k.set(idx, rows); idx = k.set(idx, cols); - idx = k.set(idx, mindisp); - idx = k.set(idx, ndisp); - idx = k.set(idx, state->preFilterCap); idx = k.set(idx, state->textureThreshold); idx = k.set(idx, state->uniquenessRatio); - idx = k.set(idx, sizeX); - idx = k.set(idx, sizeY); - idx = k.set(idx, wsz); - return k.run(3, globalThreads, localThreads, false); }