modules/gpu/src/cuda/stereobm.cu

   1 /*M///////////////////////////////////////////////////////////////////////////////////////\r
   2 //\r
   3 //  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.\r
   4 //\r
   5 //  By downloading, copying, installing or using the software you agree to this license.\r
   6 //  If you do not agree to this license, do not download, install,\r
   7 //  copy or use the software.\r
   8 //\r
   9 //\r
  10 //                           License Agreement\r
  11 //                For Open Source Computer Vision Library\r
  12 //\r
  13 // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.\r
  14 // Copyright (C) 2009, Willow Garage Inc., all rights reserved.\r
  15 // Third party copyrights are property of their respective owners.\r
  16 //\r
  17 // Redistribution and use in source and binary forms, with or without modification,\r
  18 // are permitted provided that the following conditions are met:\r
  19 //\r
  20 //   * Redistribution's of source code must retain the above copyright notice,\r
  21 //     this list of conditions and the following disclaimer.\r
  22 //\r
  23 //   * Redistribution's in binary form must reproduce the above copyright notice,\r
  24 //     this list of conditions and the following disclaimer in the documentation\r
  25 //     and/or other materials provided with the distribution.\r
  26 //\r
  27 //   * The name of the copyright holders may not be used to endorse or promote products\r
  28 //     derived from this software without specific prior written permission.\r
  29 //\r
  30 // This software is provided by the copyright holders and contributors "as is" and\r
  31 // any express or implied warranties, including, but not limited to, the implied\r
  32 // warranties of merchantability and fitness for a particular purpose are disclaimed.\r
  33 // In no event shall the Intel Corporation or contributors be liable for any direct,\r
  34 // indirect, incidental, special, exemplary, or consequential damages\r
  35 // (including, but not limited to, procurement of substitute goods or services;\r
  36 // loss of use, data, or profits; or business interruption) however caused\r
  37 // and on any theory of liability, whether in contract, strict liability,\r
  38 // or tort (including negligence or otherwise) arising in any way out of\r
  39 // the use of this software, even if advised of the possibility of such damage.\r
  40 //\r
  41 //M*/\r
  42 \r
  43 #include "internal_shared.hpp"\r
  44 \r
  45 namespace cv { namespace gpu { namespace device \r
  46 {\r
  47     namespace stereobm \r
  48     {\r
  49         //////////////////////////////////////////////////////////////////////////////////////////////////\r
  50         /////////////////////////////////////// Stereo BM ////////////////////////////////////////////////\r
  51         //////////////////////////////////////////////////////////////////////////////////////////////////\r
  52 \r
  53         #define ROWSperTHREAD 21     // the number of rows a thread will process\r
  54 \r
  55         #define BLOCK_W 128          // the thread block width (464)\r
  56         #define N_DISPARITIES 8\r
  57 \r
  58         #define STEREO_MIND 0                    // The minimum d range to check\r
  59         #define STEREO_DISP_STEP N_DISPARITIES   // the d step, must be <= 1 to avoid aliasing\r
  60 \r
  61         __constant__ unsigned int* cminSSDImage;\r
  62         __constant__ size_t cminSSD_step;\r
  63         __constant__ int cwidth;\r
  64         __constant__ int cheight;\r
  65 \r
  66         __device__ __forceinline__ int SQ(int a)\r
  67         {\r
  68             return a * a;\r
  69         }\r
  70 \r
  71         template<int RADIUS>\r
  72         __device__ unsigned int CalcSSD(volatile unsigned int *col_ssd_cache, volatile unsigned int *col_ssd)\r
  73         {       \r
  74             unsigned int cache = 0;\r
  75             unsigned int cache2 = 0;\r
  76 \r
  77             for(int i = 1; i <= RADIUS; i++)\r
  78                 cache += col_ssd[i];\r
  79 \r
  80             col_ssd_cache[0] = cache;\r
  81 \r
  82             __syncthreads();\r
  83 \r
  84             if (threadIdx.x < BLOCK_W - RADIUS)\r
  85                 cache2 = col_ssd_cache[RADIUS];\r
  86             else\r
  87                 for(int i = RADIUS + 1; i < (2 * RADIUS + 1); i++)\r
  88                     cache2 += col_ssd[i];\r
  89 \r
  90             return col_ssd[0] + cache + cache2;\r
  91         }\r
  92 \r
  93         template<int RADIUS>\r
  94         __device__ uint2 MinSSD(volatile unsigned int *col_ssd_cache, volatile unsigned int *col_ssd)\r
  95         {\r
  96             unsigned int ssd[N_DISPARITIES];\r
  97 \r
  98             //See above:  #define COL_SSD_SIZE (BLOCK_W + 2 * RADIUS)\r
  99             ssd[0] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 0 * (BLOCK_W + 2 * RADIUS));\r
 100             __syncthreads();\r
 101             ssd[1] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 1 * (BLOCK_W + 2 * RADIUS));\r
 102             __syncthreads();\r
 103             ssd[2] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 2 * (BLOCK_W + 2 * RADIUS));\r
 104             __syncthreads();\r
 105             ssd[3] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 3 * (BLOCK_W + 2 * RADIUS));\r
 106             __syncthreads();\r
 107             ssd[4] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 4 * (BLOCK_W + 2 * RADIUS));\r
 108             __syncthreads();\r
 109             ssd[5] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 5 * (BLOCK_W + 2 * RADIUS));\r
 110             __syncthreads();\r
 111             ssd[6] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 6 * (BLOCK_W + 2 * RADIUS));\r
 112             __syncthreads();\r
 113             ssd[7] = CalcSSD<RADIUS>(col_ssd_cache, col_ssd + 7 * (BLOCK_W + 2 * RADIUS));\r
 114 \r
 115             int mssd = ::min(::min(::min(ssd[0], ssd[1]), ::min(ssd[4], ssd[5])), ::min(::min(ssd[2], ssd[3]), ::min(ssd[6], ssd[7])));\r
 116 \r
 117             int bestIdx = 0;\r
 118             for (int i = 0; i < N_DISPARITIES; i++)\r
 119             {\r
 120                 if (mssd == ssd[i])\r
 121                     bestIdx = i;\r
 122             }\r
 123 \r
 124             return make_uint2(mssd, bestIdx);\r
 125         }\r
 126 \r
 127         template<int RADIUS>\r
 128         __device__ void StepDown(int idx1, int idx2, unsigned char* imageL, unsigned char* imageR, int d, volatile unsigned int *col_ssd)\r
 129         {\r
 130             unsigned char leftPixel1;\r
 131             unsigned char leftPixel2;\r
 132             unsigned char rightPixel1[8];\r
 133             unsigned char rightPixel2[8];\r
 134             unsigned int diff1, diff2;\r
 135 \r
 136             leftPixel1 = imageL[idx1];\r
 137             leftPixel2 = imageL[idx2];\r
 138 \r
 139             idx1 = idx1 - d;\r
 140             idx2 = idx2 - d;\r
 141 \r
 142             rightPixel1[7] = imageR[idx1 - 7];\r
 143             rightPixel1[0] = imageR[idx1 - 0];\r
 144             rightPixel1[1] = imageR[idx1 - 1];\r
 145             rightPixel1[2] = imageR[idx1 - 2];\r
 146             rightPixel1[3] = imageR[idx1 - 3];\r
 147             rightPixel1[4] = imageR[idx1 - 4];\r
 148             rightPixel1[5] = imageR[idx1 - 5];\r
 149             rightPixel1[6] = imageR[idx1 - 6];\r
 150 \r
 151             rightPixel2[7] = imageR[idx2 - 7];\r
 152             rightPixel2[0] = imageR[idx2 - 0];\r
 153             rightPixel2[1] = imageR[idx2 - 1];\r
 154             rightPixel2[2] = imageR[idx2 - 2];\r
 155             rightPixel2[3] = imageR[idx2 - 3];\r
 156             rightPixel2[4] = imageR[idx2 - 4];\r
 157             rightPixel2[5] = imageR[idx2 - 5];\r
 158             rightPixel2[6] = imageR[idx2 - 6];\r
 159 \r
 160             //See above:  #define COL_SSD_SIZE (BLOCK_W + 2 * RADIUS)\r
 161             diff1 = leftPixel1 - rightPixel1[0];\r
 162             diff2 = leftPixel2 - rightPixel2[0];\r
 163             col_ssd[0 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 164 \r
 165             diff1 = leftPixel1 - rightPixel1[1];\r
 166             diff2 = leftPixel2 - rightPixel2[1];\r
 167             col_ssd[1 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 168 \r
 169             diff1 = leftPixel1 - rightPixel1[2];\r
 170             diff2 = leftPixel2 - rightPixel2[2];\r
 171             col_ssd[2 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 172 \r
 173             diff1 = leftPixel1 - rightPixel1[3];\r
 174             diff2 = leftPixel2 - rightPixel2[3];\r
 175             col_ssd[3 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 176 \r
 177             diff1 = leftPixel1 - rightPixel1[4];\r
 178             diff2 = leftPixel2 - rightPixel2[4];\r
 179             col_ssd[4 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 180 \r
 181             diff1 = leftPixel1 - rightPixel1[5];\r
 182             diff2 = leftPixel2 - rightPixel2[5];\r
 183             col_ssd[5 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 184 \r
 185             diff1 = leftPixel1 - rightPixel1[6];\r
 186             diff2 = leftPixel2 - rightPixel2[6];\r
 187             col_ssd[6 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 188 \r
 189             diff1 = leftPixel1 - rightPixel1[7];\r
 190             diff2 = leftPixel2 - rightPixel2[7];\r
 191             col_ssd[7 * (BLOCK_W + 2 * RADIUS)] += SQ(diff2) - SQ(diff1);\r
 192         }\r
 193 \r
 194         template<int RADIUS>\r
 195         __device__ void InitColSSD(int x_tex, int y_tex, int im_pitch, unsigned char* imageL, unsigned char* imageR, int d, volatile unsigned int *col_ssd)\r
 196         {\r
 197             unsigned char leftPixel1;\r
 198             int idx;\r
 199             unsigned int diffa[] = {0, 0, 0, 0, 0, 0, 0, 0};\r
 200 \r
 201             for(int i = 0; i < (2 * RADIUS + 1); i++)\r
 202             {\r
 203                 idx = y_tex * im_pitch + x_tex;\r
 204                 leftPixel1 = imageL[idx];\r
 205                 idx = idx - d;\r
 206 \r
 207                 diffa[0] += SQ(leftPixel1 - imageR[idx - 0]);\r
 208                 diffa[1] += SQ(leftPixel1 - imageR[idx - 1]);\r
 209                 diffa[2] += SQ(leftPixel1 - imageR[idx - 2]);\r
 210                 diffa[3] += SQ(leftPixel1 - imageR[idx - 3]);\r
 211                 diffa[4] += SQ(leftPixel1 - imageR[idx - 4]);\r
 212                 diffa[5] += SQ(leftPixel1 - imageR[idx - 5]);\r
 213                 diffa[6] += SQ(leftPixel1 - imageR[idx - 6]);\r
 214                 diffa[7] += SQ(leftPixel1 - imageR[idx - 7]);\r
 215 \r
 216                 y_tex += 1;\r
 217             }\r
 218             //See above:  #define COL_SSD_SIZE (BLOCK_W + 2 * RADIUS)\r
 219             col_ssd[0 * (BLOCK_W + 2 * RADIUS)] = diffa[0];\r
 220             col_ssd[1 * (BLOCK_W + 2 * RADIUS)] = diffa[1];\r
 221             col_ssd[2 * (BLOCK_W + 2 * RADIUS)] = diffa[2];\r
 222             col_ssd[3 * (BLOCK_W + 2 * RADIUS)] = diffa[3];\r
 223             col_ssd[4 * (BLOCK_W + 2 * RADIUS)] = diffa[4];\r
 224             col_ssd[5 * (BLOCK_W + 2 * RADIUS)] = diffa[5];\r
 225             col_ssd[6 * (BLOCK_W + 2 * RADIUS)] = diffa[6];\r
 226             col_ssd[7 * (BLOCK_W + 2 * RADIUS)] = diffa[7];\r
 227         }\r
 228 \r
 229         template<int RADIUS>\r
 230         __global__ void stereoKernel(unsigned char *left, unsigned char *right, size_t img_step, PtrStepb disp, int maxdisp)\r
 231         {\r
 232             extern __shared__ unsigned int col_ssd_cache[];\r
 233             volatile unsigned int *col_ssd = col_ssd_cache + BLOCK_W + threadIdx.x;\r
 234             volatile unsigned int *col_ssd_extra = threadIdx.x < (2 * RADIUS) ? col_ssd + BLOCK_W : 0;  //#define N_DIRTY_PIXELS (2 * RADIUS)\r
 235 \r
 236             //#define X (blockIdx.x * BLOCK_W + threadIdx.x + STEREO_MAXD)\r
 237             int X = (blockIdx.x * BLOCK_W + threadIdx.x + maxdisp + RADIUS);\r
 238             //#define Y (__mul24(blockIdx.y, ROWSperTHREAD) + RADIUS)\r
 239             #define Y (blockIdx.y * ROWSperTHREAD + RADIUS)\r
 240             //int Y = blockIdx.y * ROWSperTHREAD + RADIUS;\r
 241 \r
 242             unsigned int* minSSDImage = cminSSDImage + X + Y * cminSSD_step;\r
 243             unsigned char* disparImage = disp.data + X + Y * disp.step;\r
 244          /*   if (X < cwidth)\r
 245             {\r
 246                 unsigned int *minSSDImage_end = minSSDImage + min(ROWSperTHREAD, cheight - Y) * minssd_step;\r
 247                 for(uint *ptr = minSSDImage; ptr != minSSDImage_end; ptr += minssd_step )\r
 248                     *ptr = 0xFFFFFFFF;\r
 249             }*/\r
 250             int end_row = ::min(ROWSperTHREAD, cheight - Y - RADIUS);\r
 251             int y_tex;\r
 252             int x_tex = X - RADIUS;\r
 253 \r
 254             if (x_tex >= cwidth)\r
 255                 return;\r
 256 \r
 257             for(int d = STEREO_MIND; d < maxdisp; d += STEREO_DISP_STEP)\r
 258             {\r
 259                 y_tex = Y - RADIUS;\r
 260 \r
 261                 InitColSSD<RADIUS>(x_tex, y_tex, img_step, left, right, d, col_ssd);\r
 262 \r
 263                 if (col_ssd_extra > 0)\r
 264                     if (x_tex + BLOCK_W < cwidth)\r
 265                         InitColSSD<RADIUS>(x_tex + BLOCK_W, y_tex, img_step, left, right, d, col_ssd_extra);\r
 266 \r
 267                 __syncthreads(); //before MinSSD function\r
 268 \r
 269                 if (X < cwidth - RADIUS && Y < cheight - RADIUS)\r
 270                 {\r
 271                     uint2 minSSD = MinSSD<RADIUS>(col_ssd_cache + threadIdx.x, col_ssd);\r
 272                     if (minSSD.x < minSSDImage[0])\r
 273                     {\r
 274                         disparImage[0] = (unsigned char)(d + minSSD.y);\r
 275                         minSSDImage[0] = minSSD.x;\r
 276                     }\r
 277                 }\r
 278 \r
 279                 for(int row = 1; row < end_row; row++)\r
 280                 {\r
 281                     int idx1 = y_tex * img_step + x_tex;\r
 282                     int idx2 = (y_tex + (2 * RADIUS + 1)) * img_step + x_tex;\r
 283 \r
 284                     __syncthreads();\r
 285 \r
 286                     StepDown<RADIUS>(idx1, idx2, left, right, d, col_ssd);\r
 287 \r
 288                     if (col_ssd_extra)\r
 289                         if (x_tex + BLOCK_W < cwidth)\r
 290                             StepDown<RADIUS>(idx1, idx2, left + BLOCK_W, right + BLOCK_W, d, col_ssd_extra);\r
 291 \r
 292                     y_tex += 1;\r
 293 \r
 294                     __syncthreads(); //before MinSSD function\r
 295 \r
 296                     if (X < cwidth - RADIUS && row < cheight - RADIUS - Y)\r
 297                     {\r
 298                         int idx = row * cminSSD_step;\r
 299                         uint2 minSSD = MinSSD<RADIUS>(col_ssd_cache + threadIdx.x, col_ssd);\r
 300                         if (minSSD.x < minSSDImage[idx])\r
 301                         {\r
 302                             disparImage[disp.step * row] = (unsigned char)(d + minSSD.y);\r
 303                             minSSDImage[idx] = minSSD.x;\r
 304                         }\r
 305                     }\r
 306                 } // for row loop\r
 307             } // for d loop\r
 308         }\r
 309 \r
 310 \r
 311         template<int RADIUS> void kernel_caller(const DevMem2Db& left, const DevMem2Db& right, const DevMem2Db& disp, int maxdisp, cudaStream_t & stream)\r
 312         {\r
 313             dim3 grid(1,1,1);\r
 314             dim3 threads(BLOCK_W, 1, 1);\r
 315 \r
 316             grid.x = divUp(left.cols - maxdisp - 2 * RADIUS, BLOCK_W);\r
 317             grid.y = divUp(left.rows - 2 * RADIUS, ROWSperTHREAD);\r
 318 \r
 319             //See above:  #define COL_SSD_SIZE (BLOCK_W + 2 * RADIUS)\r
 320             size_t smem_size = (BLOCK_W + N_DISPARITIES * (BLOCK_W + 2 * RADIUS)) * sizeof(unsigned int);\r
 321 \r
 322             stereoKernel<RADIUS><<<grid, threads, smem_size, stream>>>(left.data, right.data, left.step, disp, maxdisp);\r
 323             cudaSafeCall( cudaGetLastError() );\r
 324 \r
 325             if (stream == 0)\r
 326                 cudaSafeCall( cudaDeviceSynchronize() );\r
 327         };\r
 328 \r
 329         typedef void (*kernel_caller_t)(const DevMem2Db& left, const DevMem2Db& right, const DevMem2Db& disp, int maxdisp, cudaStream_t & stream);\r
 330 \r
 331         const static kernel_caller_t callers[] =\r
 332         {\r
 333             0,\r
 334             kernel_caller< 1>, kernel_caller< 2>, kernel_caller< 3>, kernel_caller< 4>, kernel_caller< 5>,\r
 335             kernel_caller< 6>, kernel_caller< 7>, kernel_caller< 8>, kernel_caller< 9>, kernel_caller<10>,\r
 336             kernel_caller<11>, kernel_caller<12>, kernel_caller<13>, kernel_caller<15>, kernel_caller<15>,\r
 337             kernel_caller<16>, kernel_caller<17>, kernel_caller<18>, kernel_caller<19>, kernel_caller<20>,\r
 338             kernel_caller<21>, kernel_caller<22>, kernel_caller<23>, kernel_caller<24>, kernel_caller<25>\r
 339 \r
 340             //0,0,0, 0,0,0, 0,0,kernel_caller<9>\r
 341         };\r
 342         const int calles_num = sizeof(callers)/sizeof(callers[0]);\r
 343 \r
 344         void stereoBM_GPU(const DevMem2Db& left, const DevMem2Db& right, const DevMem2Db& disp, int maxdisp, int winsz, const DevMem2D_<unsigned int>& minSSD_buf, cudaStream_t& stream)\r
 345         {\r
 346             int winsz2 = winsz >> 1;\r
 347 \r
 348             if (winsz2 == 0 || winsz2 >= calles_num)\r
 349                 cv::gpu::error("Unsupported window size", __FILE__, __LINE__, "stereoBM_GPU");\r
 350 \r
 351             //cudaSafeCall( cudaFuncSetCacheConfig(&stereoKernel, cudaFuncCachePreferL1) );\r
 352             //cudaSafeCall( cudaFuncSetCacheConfig(&stereoKernel, cudaFuncCachePreferShared) );\r
 353 \r
 354             cudaSafeCall( cudaMemset2D(disp.data, disp.step, 0, disp.cols, disp.rows) );\r
 355             cudaSafeCall( cudaMemset2D(minSSD_buf.data, minSSD_buf.step, 0xFF, minSSD_buf.cols * minSSD_buf.elemSize(), disp.rows) );\r
 356 \r
 357             cudaSafeCall( cudaMemcpyToSymbol( cwidth, &left.cols, sizeof(left.cols) ) );\r
 358             cudaSafeCall( cudaMemcpyToSymbol( cheight, &left.rows, sizeof(left.rows) ) );\r
 359             cudaSafeCall( cudaMemcpyToSymbol( cminSSDImage, &minSSD_buf.data, sizeof(minSSD_buf.data) ) );\r
 360 \r
 361             size_t minssd_step = minSSD_buf.step/minSSD_buf.elemSize();\r
 362             cudaSafeCall( cudaMemcpyToSymbol( cminSSD_step,  &minssd_step, sizeof(minssd_step) ) );\r
 363 \r
 364             callers[winsz2](left, right, disp, maxdisp, stream);\r
 365         }\r
 366 \r
 367         //////////////////////////////////////////////////////////////////////////////////////////////////\r
 368         /////////////////////////////////////// Sobel Prefiler ///////////////////////////////////////////\r
 369         //////////////////////////////////////////////////////////////////////////////////////////////////\r
 370 \r
 371         texture<unsigned char, 2, cudaReadModeElementType> texForSobel;\r
 372 \r
 373         __global__ void prefilter_kernel(DevMem2Db output, int prefilterCap)\r
 374         {\r
 375             int x = blockDim.x * blockIdx.x + threadIdx.x;\r
 376             int y = blockDim.y * blockIdx.y + threadIdx.y;\r
 377 \r
 378             if (x < output.cols && y < output.rows)\r
 379             {\r
 380                 int conv = (int)tex2D(texForSobel, x - 1, y - 1) * (-1) + (int)tex2D(texForSobel, x + 1, y - 1) * (1) +\r
 381                            (int)tex2D(texForSobel, x - 1, y    ) * (-2) + (int)tex2D(texForSobel, x + 1, y    ) * (2) +\r
 382                            (int)tex2D(texForSobel, x - 1, y + 1) * (-1) + (int)tex2D(texForSobel, x + 1, y + 1) * (1);\r
 383 \r
 384 \r
 385                 conv = ::min(::min(::max(-prefilterCap, conv), prefilterCap) + prefilterCap, 255);\r
 386                 output.ptr(y)[x] = conv & 0xFF;\r
 387             }\r
 388         }\r
 389 \r
 390         void prefilter_xsobel(const DevMem2Db& input, const DevMem2Db& output, int prefilterCap, cudaStream_t & stream)\r
 391         {\r
 392             cudaChannelFormatDesc desc = cudaCreateChannelDesc<unsigned char>();\r
 393             cudaSafeCall( cudaBindTexture2D( 0, texForSobel, input.data, desc, input.cols, input.rows, input.step ) );\r
 394 \r
 395             dim3 threads(16, 16, 1);\r
 396             dim3 grid(1, 1, 1);\r
 397 \r
 398             grid.x = divUp(input.cols, threads.x);\r
 399             grid.y = divUp(input.rows, threads.y);\r
 400 \r
 401             prefilter_kernel<<<grid, threads, 0, stream>>>(output, prefilterCap);\r
 402             cudaSafeCall( cudaGetLastError() );\r
 403 \r
 404             if (stream == 0)   \r
 405                 cudaSafeCall( cudaDeviceSynchronize() );    \r
 406 \r
 407             cudaSafeCall( cudaUnbindTexture (texForSobel ) );\r
 408         }\r
 409 \r
 410 \r
 411         //////////////////////////////////////////////////////////////////////////////////////////////////\r
 412         /////////////////////////////////// Textureness filtering ////////////////////////////////////////\r
 413         //////////////////////////////////////////////////////////////////////////////////////////////////\r
 414 \r
 415         texture<unsigned char, 2, cudaReadModeNormalizedFloat> texForTF;\r
 416 \r
 417         __device__ __forceinline__ float sobel(int x, int y)\r
 418         {\r
 419             float conv = tex2D(texForTF, x - 1, y - 1) * (-1) + tex2D(texForTF, x + 1, y - 1) * (1) +\r
 420                          tex2D(texForTF, x - 1, y    ) * (-2) + tex2D(texForTF, x + 1, y    ) * (2) +\r
 421                          tex2D(texForTF, x - 1, y + 1) * (-1) + tex2D(texForTF, x + 1, y + 1) * (1);\r
 422             return fabs(conv);\r
 423         }\r
 424 \r
 425         __device__ float CalcSums(float *cols, float *cols_cache, int winsz)\r
 426         {\r
 427             float cache = 0;\r
 428             float cache2 = 0;\r
 429             int winsz2 = winsz/2;\r
 430 \r
 431             for(int i = 1; i <= winsz2; i++)\r
 432                 cache += cols[i];\r
 433 \r
 434             cols_cache[0] = cache;\r
 435 \r
 436             __syncthreads();\r
 437 \r
 438             if (threadIdx.x < blockDim.x - winsz2)\r
 439                 cache2 = cols_cache[winsz2];\r
 440             else\r
 441                 for(int i = winsz2 + 1; i < winsz; i++)\r
 442                     cache2 += cols[i];\r
 443 \r
 444             return cols[0] + cache + cache2;\r
 445         }\r
 446 \r
 447         #define RpT (2 * ROWSperTHREAD)  // got experimentally\r
 448 \r
 449         __global__ void textureness_kernel(DevMem2Db disp, int winsz, float threshold)\r
 450         {\r
 451             int winsz2 = winsz/2;\r
 452             int n_dirty_pixels = (winsz2) * 2;\r
 453 \r
 454             extern __shared__ float cols_cache[];\r
 455             float *cols = cols_cache + blockDim.x + threadIdx.x;\r
 456             float *cols_extra = threadIdx.x < n_dirty_pixels ? cols + blockDim.x : 0;\r
 457 \r
 458             int x = blockIdx.x * blockDim.x + threadIdx.x;\r
 459             int beg_row = blockIdx.y * RpT;\r
 460             int end_row = ::min(beg_row + RpT, disp.rows);\r
 461 \r
 462             if (x < disp.cols)\r
 463             {\r
 464                 int y = beg_row;\r
 465 \r
 466                 float sum = 0;\r
 467                 float sum_extra = 0;\r
 468 \r
 469                 for(int i = y - winsz2; i <= y + winsz2; ++i)\r
 470                 {\r
 471                     sum += sobel(x - winsz2, i);\r
 472                     if (cols_extra)\r
 473                         sum_extra += sobel(x + blockDim.x - winsz2, i);\r
 474                 }\r
 475                 *cols = sum;\r
 476                 if (cols_extra)\r
 477                     *cols_extra = sum_extra;\r
 478 \r
 479                 __syncthreads();\r
 480 \r
 481                 float sum_win = CalcSums(cols, cols_cache + threadIdx.x, winsz) * 255;\r
 482                 if (sum_win < threshold)\r
 483                     disp.data[y * disp.step + x] = 0;\r
 484 \r
 485                 __syncthreads();\r
 486 \r
 487                 for(int y = beg_row + 1; y < end_row; ++y)\r
 488                 {\r
 489                     sum = sum - sobel(x - winsz2, y - winsz2 - 1) + sobel(x - winsz2, y + winsz2);\r
 490                     *cols = sum;\r
 491 \r
 492                     if (cols_extra)\r
 493                     {\r
 494                         sum_extra = sum_extra - sobel(x + blockDim.x - winsz2, y - winsz2 - 1) + sobel(x + blockDim.x - winsz2, y + winsz2);\r
 495                         *cols_extra = sum_extra;\r
 496                     }\r
 497 \r
 498                     __syncthreads();\r
 499                     float sum_win = CalcSums(cols, cols_cache + threadIdx.x, winsz) * 255;\r
 500                     if (sum_win < threshold)\r
 501                         disp.data[y * disp.step + x] = 0;\r
 502 \r
 503                     __syncthreads();\r
 504                 }\r
 505             }\r
 506         }\r
 507 \r
 508         void postfilter_textureness(const DevMem2Db& input, int winsz, float avgTexturenessThreshold, const DevMem2Db& disp, cudaStream_t & stream)\r
 509         {\r
 510             avgTexturenessThreshold *= winsz * winsz;\r
 511 \r
 512             texForTF.filterMode     = cudaFilterModeLinear;\r
 513             texForTF.addressMode[0] = cudaAddressModeWrap;\r
 514             texForTF.addressMode[1] = cudaAddressModeWrap;\r
 515 \r
 516             cudaChannelFormatDesc desc = cudaCreateChannelDesc<unsigned char>();\r
 517             cudaSafeCall( cudaBindTexture2D( 0, texForTF, input.data, desc, input.cols, input.rows, input.step ) );\r
 518 \r
 519             dim3 threads(128, 1, 1);\r
 520             dim3 grid(1, 1, 1);\r
 521 \r
 522             grid.x = divUp(input.cols, threads.x);\r
 523             grid.y = divUp(input.rows, RpT);\r
 524 \r
 525             size_t smem_size = (threads.x + threads.x + (winsz/2) * 2 ) * sizeof(float);\r
 526             textureness_kernel<<<grid, threads, smem_size, stream>>>(disp, winsz, avgTexturenessThreshold);\r
 527             cudaSafeCall( cudaGetLastError() );\r
 528 \r
 529             if (stream == 0)\r
 530                 cudaSafeCall( cudaDeviceSynchronize() );\r
 531 \r
 532             cudaSafeCall( cudaUnbindTexture (texForTF) );\r
 533         }\r
 534     } // namespace stereobm\r
 535 }}} // namespace cv { namespace gpu { namespace device\r