modules/gpu/src/cuda/matrix_reductions.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 #include "opencv2/gpu/device/limits.hpp"\r
  45 #include "opencv2/gpu/device/saturate_cast.hpp"\r
  46 #include "opencv2/gpu/device/vec_math.hpp"\r
  47 \r
  48 namespace cv { namespace gpu { namespace device \r
  49 {\r
  50     namespace matrix_reductions \r
  51     {\r
  52         // Performs reduction in shared memory\r
  53         template <int size, typename T>\r
  54         __device__ void sumInSmem(volatile T* data, const uint tid)\r
  55         {\r
  56             T sum = data[tid];\r
  57 \r
  58             if (size >= 512) { if (tid < 256) { data[tid] = sum = sum + data[tid + 256]; } __syncthreads(); }\r
  59             if (size >= 256) { if (tid < 128) { data[tid] = sum = sum + data[tid + 128]; } __syncthreads(); }\r
  60             if (size >= 128) { if (tid < 64) { data[tid] = sum = sum + data[tid + 64]; } __syncthreads(); }\r
  61 \r
  62             if (tid < 32)\r
  63             {\r
  64                 if (size >= 64) data[tid] = sum = sum + data[tid + 32];\r
  65                 if (size >= 32) data[tid] = sum = sum + data[tid + 16];\r
  66                 if (size >= 16) data[tid] = sum = sum + data[tid + 8];\r
  67                 if (size >= 8) data[tid] = sum = sum + data[tid + 4];\r
  68                 if (size >= 4) data[tid] = sum = sum + data[tid + 2];\r
  69                 if (size >= 2) data[tid] = sum = sum + data[tid + 1];\r
  70             }\r
  71         }\r
  72 \r
  73         struct Mask8U\r
  74         {\r
  75             explicit Mask8U(PtrStepb mask): mask(mask) {}\r
  76 \r
  77             __device__ __forceinline__ bool operator()(int y, int x) const \r
  78             { \r
  79                 return mask.ptr(y)[x]; \r
  80             }\r
  81 \r
  82             PtrStepb mask;\r
  83         };\r
  84 \r
  85         struct MaskTrue \r
  86         { \r
  87             __device__ __forceinline__ bool operator()(int y, int x) const \r
  88             { \r
  89                 return true; \r
  90             }\r
  91             __device__ __forceinline__ MaskTrue(){}\r
  92             __device__ __forceinline__ MaskTrue(const MaskTrue& mask_){}\r
  93         };\r
  94 \r
  95         //////////////////////////////////////////////////////////////////////////////\r
  96         // Min max\r
  97 \r
  98         // To avoid shared bank conflicts we convert each value into value of \r
  99         // appropriate type (32 bits minimum)\r
 100         template <typename T> struct MinMaxTypeTraits {};\r
 101         template <> struct MinMaxTypeTraits<uchar> { typedef int best_type; };\r
 102         template <> struct MinMaxTypeTraits<char> { typedef int best_type; };\r
 103         template <> struct MinMaxTypeTraits<ushort> { typedef int best_type; };\r
 104         template <> struct MinMaxTypeTraits<short> { typedef int best_type; };\r
 105         template <> struct MinMaxTypeTraits<int> { typedef int best_type; };\r
 106         template <> struct MinMaxTypeTraits<float> { typedef float best_type; };\r
 107         template <> struct MinMaxTypeTraits<double> { typedef double best_type; };\r
 108 \r
 109         namespace minmax \r
 110         {\r
 111             __constant__ int ctwidth;\r
 112             __constant__ int ctheight;\r
 113 \r
 114             // Global counter of blocks finished its work\r
 115             __device__ uint blocks_finished = 0;\r
 116 \r
 117 \r
 118             // Estimates good thread configuration\r
 119             //  - threads variable satisfies to threads.x * threads.y == 256\r
 120             void estimateThreadCfg(int cols, int rows, dim3& threads, dim3& grid)\r
 121             {\r
 122                 threads = dim3(32, 8);\r
 123                 grid = dim3(divUp(cols, threads.x * 8), divUp(rows, threads.y * 32));\r
 124                 grid.x = std::min(grid.x, threads.x);\r
 125                 grid.y = std::min(grid.y, threads.y);\r
 126             }\r
 127 \r
 128 \r
 129             // Returns required buffer sizes\r
 130             void getBufSizeRequired(int cols, int rows, int elem_size, int& bufcols, int& bufrows)\r
 131             {\r
 132                 dim3 threads, grid;\r
 133                 estimateThreadCfg(cols, rows, threads, grid);\r
 134                 bufcols = grid.x * grid.y * elem_size; \r
 135                 bufrows = 2;\r
 136             }\r
 137 \r
 138 \r
 139             // Estimates device constants which are used in the kernels using specified thread configuration\r
 140             void setKernelConsts(int cols, int rows, const dim3& threads, const dim3& grid)\r
 141             {        \r
 142                 int twidth = divUp(divUp(cols, grid.x), threads.x);\r
 143                 int theight = divUp(divUp(rows, grid.y), threads.y);\r
 144                 cudaSafeCall(cudaMemcpyToSymbol(ctwidth, &twidth, sizeof(ctwidth))); \r
 145                 cudaSafeCall(cudaMemcpyToSymbol(ctheight, &theight, sizeof(ctheight))); \r
 146             }  \r
 147 \r
 148 \r
 149             // Does min and max in shared memory\r
 150             template <typename T>\r
 151             __device__ __forceinline__ void merge(uint tid, uint offset, volatile T* minval, volatile T* maxval)\r
 152             {\r
 153                 minval[tid] = ::min(minval[tid], minval[tid + offset]);\r
 154                 maxval[tid] = ::max(maxval[tid], maxval[tid + offset]);\r
 155             }\r
 156 \r
 157 \r
 158             template <int size, typename T>\r
 159             __device__ void findMinMaxInSmem(volatile T* minval, volatile T* maxval, const uint tid)\r
 160             {\r
 161                 if (size >= 512) { if (tid < 256) { merge(tid, 256, minval, maxval); } __syncthreads(); }\r
 162                 if (size >= 256) { if (tid < 128) { merge(tid, 128, minval, maxval); }  __syncthreads(); }\r
 163                 if (size >= 128) { if (tid < 64) { merge(tid, 64, minval, maxval); } __syncthreads(); }\r
 164 \r
 165                 if (tid < 32)\r
 166                 {\r
 167                     if (size >= 64) merge(tid, 32, minval, maxval);\r
 168                     if (size >= 32) merge(tid, 16, minval, maxval);\r
 169                     if (size >= 16) merge(tid, 8, minval, maxval);\r
 170                     if (size >= 8) merge(tid, 4, minval, maxval);\r
 171                     if (size >= 4) merge(tid, 2, minval, maxval);\r
 172                     if (size >= 2) merge(tid, 1, minval, maxval);\r
 173                 }\r
 174             }\r
 175 \r
 176 \r
 177             template <int nthreads, typename T, typename Mask>\r
 178             __global__ void minMaxKernel(const DevMem2Db src, Mask mask, T* minval, T* maxval)\r
 179             {\r
 180                 typedef typename MinMaxTypeTraits<T>::best_type best_type;\r
 181                 __shared__ best_type sminval[nthreads];\r
 182                 __shared__ best_type smaxval[nthreads];\r
 183 \r
 184                 uint x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
 185                 uint y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
 186                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 187 \r
 188                 T mymin = numeric_limits<T>::max();\r
 189                 T mymax = numeric_limits<T>::is_signed ? -numeric_limits<T>::max() : numeric_limits<T>::min();\r
 190                 uint y_end = ::min(y0 + (ctheight - 1) * blockDim.y + 1, src.rows);\r
 191                 uint x_end = ::min(x0 + (ctwidth - 1) * blockDim.x + 1, src.cols);\r
 192                 for (uint y = y0; y < y_end; y += blockDim.y)\r
 193                 {\r
 194                     const T* src_row = (const T*)src.ptr(y);\r
 195                     for (uint x = x0; x < x_end; x += blockDim.x)\r
 196                     {\r
 197                         T val = src_row[x];\r
 198                         if (mask(y, x)) \r
 199                         { \r
 200                             mymin = ::min(mymin, val); \r
 201                             mymax = ::max(mymax, val); \r
 202                         }\r
 203                     }\r
 204                 }\r
 205 \r
 206                 sminval[tid] = mymin;\r
 207                 smaxval[tid] = mymax;\r
 208                 __syncthreads();\r
 209 \r
 210                 findMinMaxInSmem<nthreads, best_type>(sminval, smaxval, tid);\r
 211 \r
 212                 if (tid == 0) \r
 213                 {\r
 214                     minval[blockIdx.y * gridDim.x + blockIdx.x] = (T)sminval[0];\r
 215                     maxval[blockIdx.y * gridDim.x + blockIdx.x] = (T)smaxval[0];\r
 216                 }\r
 217 \r
 218             #if __CUDA_ARCH__ >= 110\r
 219                         __shared__ bool is_last;\r
 220 \r
 221                         if (tid == 0)\r
 222                         {\r
 223                                 minval[blockIdx.y * gridDim.x + blockIdx.x] = (T)sminval[0];\r
 224                     maxval[blockIdx.y * gridDim.x + blockIdx.x] = (T)smaxval[0];\r
 225                                 __threadfence();\r
 226 \r
 227                                 uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
 228                                 is_last = ticket == gridDim.x * gridDim.y - 1;\r
 229                         }\r
 230 \r
 231                         __syncthreads();\r
 232 \r
 233                         if (is_last)\r
 234                         {\r
 235                     uint idx = ::min(tid, gridDim.x * gridDim.y - 1);\r
 236 \r
 237                     sminval[tid] = minval[idx];\r
 238                     smaxval[tid] = maxval[idx];\r
 239                     __syncthreads();\r
 240 \r
 241                                 findMinMaxInSmem<nthreads, best_type>(sminval, smaxval, tid);\r
 242 \r
 243                     if (tid == 0) \r
 244                     {\r
 245                         minval[0] = (T)sminval[0];\r
 246                         maxval[0] = (T)smaxval[0];\r
 247                         blocks_finished = 0;\r
 248                     }\r
 249                         }\r
 250             #else\r
 251                 if (tid == 0) \r
 252                 {\r
 253                     minval[blockIdx.y * gridDim.x + blockIdx.x] = (T)sminval[0];\r
 254                     maxval[blockIdx.y * gridDim.x + blockIdx.x] = (T)smaxval[0];\r
 255                 }\r
 256             #endif\r
 257             }\r
 258 \r
 259    \r
 260             template <typename T>\r
 261             void minMaxMaskCaller(const DevMem2Db src, const PtrStepb mask, double* minval, double* maxval, PtrStepb buf)\r
 262             {\r
 263                 dim3 threads, grid;\r
 264                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 265                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 266 \r
 267                 T* minval_buf = (T*)buf.ptr(0);\r
 268                 T* maxval_buf = (T*)buf.ptr(1);\r
 269 \r
 270                 minMaxKernel<256, T, Mask8U><<<grid, threads>>>(src, Mask8U(mask), minval_buf, maxval_buf);\r
 271                 cudaSafeCall( cudaGetLastError() );\r
 272 \r
 273                 cudaSafeCall( cudaDeviceSynchronize() );\r
 274 \r
 275                 T minval_, maxval_;\r
 276                 cudaSafeCall( cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 277                 cudaSafeCall( cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 278                 *minval = minval_;\r
 279                 *maxval = maxval_;\r
 280             }  \r
 281 \r
 282             template void minMaxMaskCaller<uchar>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 283             template void minMaxMaskCaller<char>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 284             template void minMaxMaskCaller<ushort>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 285             template void minMaxMaskCaller<short>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 286             template void minMaxMaskCaller<int>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 287             template void minMaxMaskCaller<float>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 288             template void minMaxMaskCaller<double>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 289 \r
 290 \r
 291             template <typename T>\r
 292             void minMaxCaller(const DevMem2Db src, double* minval, double* maxval, PtrStepb buf)\r
 293             {\r
 294                 dim3 threads, grid;\r
 295                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 296                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 297 \r
 298                 T* minval_buf = (T*)buf.ptr(0);\r
 299                 T* maxval_buf = (T*)buf.ptr(1);\r
 300 \r
 301                 minMaxKernel<256, T, MaskTrue><<<grid, threads>>>(src, MaskTrue(), minval_buf, maxval_buf);\r
 302                 cudaSafeCall( cudaGetLastError() );\r
 303 \r
 304                 cudaSafeCall( cudaDeviceSynchronize() );\r
 305 \r
 306                 T minval_, maxval_;\r
 307                 cudaSafeCall( cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 308                 cudaSafeCall( cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 309                 *minval = minval_;\r
 310                 *maxval = maxval_;\r
 311             }  \r
 312 \r
 313             template void minMaxCaller<uchar>(const DevMem2Db, double*, double*, PtrStepb);\r
 314             template void minMaxCaller<char>(const DevMem2Db, double*, double*, PtrStepb);\r
 315             template void minMaxCaller<ushort>(const DevMem2Db, double*, double*, PtrStepb);\r
 316             template void minMaxCaller<short>(const DevMem2Db, double*, double*, PtrStepb);\r
 317             template void minMaxCaller<int>(const DevMem2Db, double*, double*, PtrStepb);\r
 318             template void minMaxCaller<float>(const DevMem2Db, double*,double*, PtrStepb);\r
 319             template void minMaxCaller<double>(const DevMem2Db, double*, double*, PtrStepb);\r
 320 \r
 321 \r
 322             template <int nthreads, typename T>\r
 323             __global__ void minMaxPass2Kernel(T* minval, T* maxval, int size)\r
 324             {\r
 325                 typedef typename MinMaxTypeTraits<T>::best_type best_type;\r
 326                 __shared__ best_type sminval[nthreads];\r
 327                 __shared__ best_type smaxval[nthreads];\r
 328                 \r
 329                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 330                 uint idx = ::min(tid, size - 1);\r
 331 \r
 332                 sminval[tid] = minval[idx];\r
 333                 smaxval[tid] = maxval[idx];\r
 334                 __syncthreads();\r
 335 \r
 336                 findMinMaxInSmem<nthreads, best_type>(sminval, smaxval, tid);\r
 337 \r
 338                 if (tid == 0) \r
 339                 {\r
 340                     minval[0] = (T)sminval[0];\r
 341                     maxval[0] = (T)smaxval[0];\r
 342                 }\r
 343             }\r
 344 \r
 345 \r
 346             template <typename T>\r
 347             void minMaxMaskMultipassCaller(const DevMem2Db src, const PtrStepb mask, double* minval, double* maxval, PtrStepb buf)\r
 348             {\r
 349                 dim3 threads, grid;\r
 350                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 351                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 352 \r
 353                 T* minval_buf = (T*)buf.ptr(0);\r
 354                 T* maxval_buf = (T*)buf.ptr(1);\r
 355 \r
 356                 minMaxKernel<256, T, Mask8U><<<grid, threads>>>(src, Mask8U(mask), minval_buf, maxval_buf);\r
 357                 cudaSafeCall( cudaGetLastError() );\r
 358                 minMaxPass2Kernel<256, T><<<1, 256>>>(minval_buf, maxval_buf, grid.x * grid.y);\r
 359                 cudaSafeCall( cudaGetLastError() );\r
 360 \r
 361                 cudaSafeCall(cudaDeviceSynchronize());\r
 362 \r
 363                 T minval_, maxval_;\r
 364                 cudaSafeCall( cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 365                 cudaSafeCall( cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 366                 *minval = minval_;\r
 367                 *maxval = maxval_;\r
 368             }\r
 369 \r
 370             template void minMaxMaskMultipassCaller<uchar>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 371             template void minMaxMaskMultipassCaller<char>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 372             template void minMaxMaskMultipassCaller<ushort>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 373             template void minMaxMaskMultipassCaller<short>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 374             template void minMaxMaskMultipassCaller<int>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 375             template void minMaxMaskMultipassCaller<float>(const DevMem2Db, const PtrStepb, double*, double*, PtrStepb);\r
 376 \r
 377 \r
 378             template <typename T>\r
 379             void minMaxMultipassCaller(const DevMem2Db src, double* minval, double* maxval, PtrStepb buf)\r
 380             {\r
 381                 dim3 threads, grid;\r
 382                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 383                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 384 \r
 385                 T* minval_buf = (T*)buf.ptr(0);\r
 386                 T* maxval_buf = (T*)buf.ptr(1);\r
 387 \r
 388                 minMaxKernel<256, T, MaskTrue><<<grid, threads>>>(src, MaskTrue(), minval_buf, maxval_buf);\r
 389                 cudaSafeCall( cudaGetLastError() );\r
 390                 minMaxPass2Kernel<256, T><<<1, 256>>>(minval_buf, maxval_buf, grid.x * grid.y);\r
 391                 cudaSafeCall( cudaGetLastError() );\r
 392 \r
 393                 cudaSafeCall( cudaDeviceSynchronize() );\r
 394 \r
 395                 T minval_, maxval_;\r
 396                 cudaSafeCall( cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 397                 cudaSafeCall( cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 398                 *minval = minval_;\r
 399                 *maxval = maxval_;\r
 400             }\r
 401 \r
 402             template void minMaxMultipassCaller<uchar>(const DevMem2Db, double*, double*, PtrStepb);\r
 403             template void minMaxMultipassCaller<char>(const DevMem2Db, double*, double*, PtrStepb);\r
 404             template void minMaxMultipassCaller<ushort>(const DevMem2Db, double*, double*, PtrStepb);\r
 405             template void minMaxMultipassCaller<short>(const DevMem2Db, double*, double*, PtrStepb);\r
 406             template void minMaxMultipassCaller<int>(const DevMem2Db, double*, double*, PtrStepb);\r
 407             template void minMaxMultipassCaller<float>(const DevMem2Db, double*, double*, PtrStepb);\r
 408         } // namespace minmax\r
 409 \r
 410         ///////////////////////////////////////////////////////////////////////////////\r
 411         // minMaxLoc\r
 412 \r
 413         namespace minmaxloc \r
 414         {\r
 415             __constant__ int ctwidth;\r
 416             __constant__ int ctheight;\r
 417 \r
 418             // Global counter of blocks finished its work\r
 419             __device__ uint blocks_finished = 0;\r
 420 \r
 421 \r
 422             // Estimates good thread configuration\r
 423             //  - threads variable satisfies to threads.x * threads.y == 256\r
 424             void estimateThreadCfg(int cols, int rows, dim3& threads, dim3& grid)\r
 425             {\r
 426                 threads = dim3(32, 8);\r
 427                 grid = dim3(divUp(cols, threads.x * 8), divUp(rows, threads.y * 32));\r
 428                 grid.x = std::min(grid.x, threads.x);\r
 429                 grid.y = std::min(grid.y, threads.y);\r
 430             }\r
 431 \r
 432 \r
 433             // Returns required buffer sizes\r
 434             void getBufSizeRequired(int cols, int rows, int elem_size, int& b1cols, \r
 435                                     int& b1rows, int& b2cols, int& b2rows)\r
 436             {\r
 437                 dim3 threads, grid;\r
 438                 estimateThreadCfg(cols, rows, threads, grid);\r
 439                 b1cols = grid.x * grid.y * elem_size; // For values\r
 440                 b1rows = 2;\r
 441                 b2cols = grid.x * grid.y * sizeof(int); // For locations\r
 442                 b2rows = 2;\r
 443             }\r
 444 \r
 445 \r
 446             // Estimates device constants which are used in the kernels using specified thread configuration\r
 447             void setKernelConsts(int cols, int rows, const dim3& threads, const dim3& grid)\r
 448             {        \r
 449                 int twidth = divUp(divUp(cols, grid.x), threads.x);\r
 450                 int theight = divUp(divUp(rows, grid.y), threads.y);\r
 451                 cudaSafeCall(cudaMemcpyToSymbol(ctwidth, &twidth, sizeof(ctwidth))); \r
 452                 cudaSafeCall(cudaMemcpyToSymbol(ctheight, &theight, sizeof(ctheight))); \r
 453             }  \r
 454 \r
 455 \r
 456             template <typename T>\r
 457             __device__ void merge(uint tid, uint offset, volatile T* minval, volatile T* maxval, \r
 458                                   volatile uint* minloc, volatile uint* maxloc)\r
 459             {\r
 460                 T val = minval[tid + offset];\r
 461                 if (val < minval[tid])\r
 462                 {\r
 463                     minval[tid] = val;\r
 464                     minloc[tid] = minloc[tid + offset];\r
 465                 }\r
 466                 val = maxval[tid + offset];\r
 467                 if (val > maxval[tid])\r
 468                 {\r
 469                     maxval[tid] = val;\r
 470                     maxloc[tid] = maxloc[tid + offset];\r
 471                 }\r
 472             }\r
 473 \r
 474 \r
 475             template <int size, typename T>\r
 476             __device__ void findMinMaxLocInSmem(volatile T* minval, volatile T* maxval, volatile uint* minloc, \r
 477                                                 volatile uint* maxloc, const uint tid)\r
 478             {\r
 479                 if (size >= 512) { if (tid < 256) { merge(tid, 256, minval, maxval, minloc, maxloc); } __syncthreads(); }\r
 480                 if (size >= 256) { if (tid < 128) { merge(tid, 128, minval, maxval, minloc, maxloc); }  __syncthreads(); }\r
 481                 if (size >= 128) { if (tid < 64) { merge(tid, 64, minval, maxval, minloc, maxloc); } __syncthreads(); }\r
 482 \r
 483                 if (tid < 32)\r
 484                 {\r
 485                     if (size >= 64) merge(tid, 32, minval, maxval, minloc, maxloc);\r
 486                     if (size >= 32) merge(tid, 16, minval, maxval, minloc, maxloc);\r
 487                     if (size >= 16) merge(tid, 8, minval, maxval, minloc, maxloc);\r
 488                     if (size >= 8) merge(tid, 4, minval, maxval, minloc, maxloc);\r
 489                     if (size >= 4) merge(tid, 2, minval, maxval, minloc, maxloc);\r
 490                     if (size >= 2) merge(tid, 1, minval, maxval, minloc, maxloc);\r
 491                 }\r
 492             }\r
 493 \r
 494 \r
 495             template <int nthreads, typename T, typename Mask>\r
 496             __global__ void minMaxLocKernel(const DevMem2Db src, Mask mask, T* minval, T* maxval, \r
 497                                             uint* minloc, uint* maxloc)\r
 498             {\r
 499                 typedef typename MinMaxTypeTraits<T>::best_type best_type;\r
 500                 __shared__ best_type sminval[nthreads];\r
 501                 __shared__ best_type smaxval[nthreads];\r
 502                 __shared__ uint sminloc[nthreads];\r
 503                 __shared__ uint smaxloc[nthreads];\r
 504 \r
 505                 uint x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
 506                 uint y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
 507                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 508 \r
 509                 T mymin = numeric_limits<T>::max();\r
 510                 T mymax = numeric_limits<T>::is_signed ? -numeric_limits<T>::max() : numeric_limits<T>::min(); \r
 511                 uint myminloc = 0;\r
 512                 uint mymaxloc = 0;\r
 513                 uint y_end = ::min(y0 + (ctheight - 1) * blockDim.y + 1, src.rows);\r
 514                 uint x_end = ::min(x0 + (ctwidth - 1) * blockDim.x + 1, src.cols);\r
 515 \r
 516                 for (uint y = y0; y < y_end; y += blockDim.y)\r
 517                 {\r
 518                     const T* ptr = (const T*)src.ptr(y);\r
 519                     for (uint x = x0; x < x_end; x += blockDim.x)\r
 520                     {\r
 521                         if (mask(y, x))\r
 522                         {\r
 523                             T val = ptr[x];\r
 524                             if (val <= mymin) { mymin = val; myminloc = y * src.cols + x; }\r
 525                             if (val >= mymax) { mymax = val; mymaxloc = y * src.cols + x; }\r
 526                         }\r
 527                     }\r
 528                 }\r
 529 \r
 530                 sminval[tid] = mymin; \r
 531                 smaxval[tid] = mymax;\r
 532                 sminloc[tid] = myminloc;\r
 533                 smaxloc[tid] = mymaxloc;\r
 534                 __syncthreads();\r
 535 \r
 536                 findMinMaxLocInSmem<nthreads, best_type>(sminval, smaxval, sminloc, smaxloc, tid);\r
 537 \r
 538             #if __CUDA_ARCH__ >= 110\r
 539                         __shared__ bool is_last;\r
 540 \r
 541                         if (tid == 0)\r
 542                         {\r
 543                                 minval[blockIdx.y * gridDim.x + blockIdx.x] = (T)sminval[0];\r
 544                     maxval[blockIdx.y * gridDim.x + blockIdx.x] = (T)smaxval[0];\r
 545                     minloc[blockIdx.y * gridDim.x + blockIdx.x] = sminloc[0];\r
 546                     maxloc[blockIdx.y * gridDim.x + blockIdx.x] = smaxloc[0];\r
 547                                 __threadfence();\r
 548 \r
 549                                 uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
 550                                 is_last = ticket == gridDim.x * gridDim.y - 1;\r
 551                         }\r
 552 \r
 553                         __syncthreads();\r
 554 \r
 555                         if (is_last)\r
 556                         {\r
 557                     uint idx = ::min(tid, gridDim.x * gridDim.y - 1);\r
 558 \r
 559                     sminval[tid] = minval[idx];\r
 560                     smaxval[tid] = maxval[idx];\r
 561                     sminloc[tid] = minloc[idx];\r
 562                     smaxloc[tid] = maxloc[idx];\r
 563                     __syncthreads();\r
 564 \r
 565                                 findMinMaxLocInSmem<nthreads, best_type>(sminval, smaxval, sminloc, smaxloc, tid);\r
 566 \r
 567                     if (tid == 0) \r
 568                     {\r
 569                         minval[0] = (T)sminval[0];\r
 570                         maxval[0] = (T)smaxval[0];\r
 571                         minloc[0] = sminloc[0];\r
 572                         maxloc[0] = smaxloc[0];\r
 573                         blocks_finished = 0;\r
 574                     }\r
 575                         }\r
 576             #else\r
 577                 if (tid == 0) \r
 578                 {\r
 579                     minval[blockIdx.y * gridDim.x + blockIdx.x] = (T)sminval[0];\r
 580                     maxval[blockIdx.y * gridDim.x + blockIdx.x] = (T)smaxval[0];\r
 581                     minloc[blockIdx.y * gridDim.x + blockIdx.x] = sminloc[0];\r
 582                     maxloc[blockIdx.y * gridDim.x + blockIdx.x] = smaxloc[0];\r
 583                 }\r
 584             #endif\r
 585             }\r
 586 \r
 587 \r
 588             template <typename T>\r
 589             void minMaxLocMaskCaller(const DevMem2Db src, const PtrStepb mask, double* minval, double* maxval, \r
 590                                      int minloc[2], int maxloc[2], PtrStepb valbuf, PtrStepb locbuf)\r
 591             {\r
 592                 dim3 threads, grid;\r
 593                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 594                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 595 \r
 596                 T* minval_buf = (T*)valbuf.ptr(0);\r
 597                 T* maxval_buf = (T*)valbuf.ptr(1);\r
 598                 uint* minloc_buf = (uint*)locbuf.ptr(0);\r
 599                 uint* maxloc_buf = (uint*)locbuf.ptr(1);\r
 600 \r
 601                 minMaxLocKernel<256, T, Mask8U><<<grid, threads>>>(src, Mask8U(mask), minval_buf, maxval_buf, \r
 602                                                                    minloc_buf, maxloc_buf);\r
 603                 cudaSafeCall( cudaGetLastError() );\r
 604 \r
 605                 cudaSafeCall( cudaDeviceSynchronize() );\r
 606 \r
 607                 T minval_, maxval_;\r
 608                 cudaSafeCall( cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 609                 cudaSafeCall( cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost) );\r
 610                 *minval = minval_;\r
 611                 *maxval = maxval_;\r
 612 \r
 613                 uint minloc_, maxloc_;\r
 614                 cudaSafeCall( cudaMemcpy(&minloc_, minloc_buf, sizeof(int), cudaMemcpyDeviceToHost) );\r
 615                 cudaSafeCall( cudaMemcpy(&maxloc_, maxloc_buf, sizeof(int), cudaMemcpyDeviceToHost) );\r
 616                 minloc[1] = minloc_ / src.cols; minloc[0] = minloc_ - minloc[1] * src.cols;\r
 617                 maxloc[1] = maxloc_ / src.cols; maxloc[0] = maxloc_ - maxloc[1] * src.cols;\r
 618             }\r
 619 \r
 620             template void minMaxLocMaskCaller<uchar>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 621             template void minMaxLocMaskCaller<char>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 622             template void minMaxLocMaskCaller<ushort>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 623             template void minMaxLocMaskCaller<short>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 624             template void minMaxLocMaskCaller<int>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 625             template void minMaxLocMaskCaller<float>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 626             template void minMaxLocMaskCaller<double>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 627 \r
 628 \r
 629             template <typename T>\r
 630             void minMaxLocCaller(const DevMem2Db src, double* minval, double* maxval, \r
 631                                  int minloc[2], int maxloc[2], PtrStepb valbuf, PtrStepb locbuf)\r
 632             {\r
 633                 dim3 threads, grid;\r
 634                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 635                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 636 \r
 637                 T* minval_buf = (T*)valbuf.ptr(0);\r
 638                 T* maxval_buf = (T*)valbuf.ptr(1);\r
 639                 uint* minloc_buf = (uint*)locbuf.ptr(0);\r
 640                 uint* maxloc_buf = (uint*)locbuf.ptr(1);\r
 641 \r
 642                 minMaxLocKernel<256, T, MaskTrue><<<grid, threads>>>(src, MaskTrue(), minval_buf, maxval_buf, \r
 643                                                                      minloc_buf, maxloc_buf);\r
 644                 cudaSafeCall( cudaGetLastError() );\r
 645 \r
 646                 cudaSafeCall( cudaDeviceSynchronize() );\r
 647 \r
 648                 T minval_, maxval_;\r
 649                 cudaSafeCall(cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 650                 cudaSafeCall(cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 651                 *minval = minval_;\r
 652                 *maxval = maxval_;\r
 653 \r
 654                 uint minloc_, maxloc_;\r
 655                 cudaSafeCall(cudaMemcpy(&minloc_, minloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 656                 cudaSafeCall(cudaMemcpy(&maxloc_, maxloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 657                 minloc[1] = minloc_ / src.cols; minloc[0] = minloc_ - minloc[1] * src.cols;\r
 658                 maxloc[1] = maxloc_ / src.cols; maxloc[0] = maxloc_ - maxloc[1] * src.cols;\r
 659             }\r
 660 \r
 661             template void minMaxLocCaller<uchar>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 662             template void minMaxLocCaller<char>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 663             template void minMaxLocCaller<ushort>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 664             template void minMaxLocCaller<short>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 665             template void minMaxLocCaller<int>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 666             template void minMaxLocCaller<float>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 667             template void minMaxLocCaller<double>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 668 \r
 669 \r
 670             // This kernel will be used only when compute capability is 1.0\r
 671             template <int nthreads, typename T>\r
 672             __global__ void minMaxLocPass2Kernel(T* minval, T* maxval, uint* minloc, uint* maxloc, int size)\r
 673             {\r
 674                 typedef typename MinMaxTypeTraits<T>::best_type best_type;\r
 675                 __shared__ best_type sminval[nthreads];\r
 676                 __shared__ best_type smaxval[nthreads];\r
 677                 __shared__ uint sminloc[nthreads];\r
 678                 __shared__ uint smaxloc[nthreads];\r
 679 \r
 680                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 681                 uint idx = ::min(tid, size - 1);\r
 682 \r
 683                 sminval[tid] = minval[idx];\r
 684                 smaxval[tid] = maxval[idx];\r
 685                 sminloc[tid] = minloc[idx];\r
 686                 smaxloc[tid] = maxloc[idx];\r
 687                 __syncthreads();\r
 688 \r
 689                 findMinMaxLocInSmem<nthreads, best_type>(sminval, smaxval, sminloc, smaxloc, tid);\r
 690 \r
 691                 if (tid == 0) \r
 692                 {\r
 693                     minval[0] = (T)sminval[0];\r
 694                     maxval[0] = (T)smaxval[0];\r
 695                     minloc[0] = sminloc[0];\r
 696                     maxloc[0] = smaxloc[0];\r
 697                 }\r
 698             }\r
 699 \r
 700 \r
 701             template <typename T>\r
 702             void minMaxLocMaskMultipassCaller(const DevMem2Db src, const PtrStepb mask, double* minval, double* maxval, \r
 703                                               int minloc[2], int maxloc[2], PtrStepb valbuf, PtrStepb locbuf)\r
 704             {\r
 705                 dim3 threads, grid;\r
 706                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 707                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 708 \r
 709                 T* minval_buf = (T*)valbuf.ptr(0);\r
 710                 T* maxval_buf = (T*)valbuf.ptr(1);\r
 711                 uint* minloc_buf = (uint*)locbuf.ptr(0);\r
 712                 uint* maxloc_buf = (uint*)locbuf.ptr(1);\r
 713 \r
 714                 minMaxLocKernel<256, T, Mask8U><<<grid, threads>>>(src, Mask8U(mask), minval_buf, maxval_buf, \r
 715                                                                    minloc_buf, maxloc_buf);\r
 716                 cudaSafeCall( cudaGetLastError() );\r
 717                 minMaxLocPass2Kernel<256, T><<<1, 256>>>(minval_buf, maxval_buf, minloc_buf, maxloc_buf, grid.x * grid.y);\r
 718                 cudaSafeCall( cudaGetLastError() );\r
 719 \r
 720                 cudaSafeCall( cudaDeviceSynchronize() );\r
 721 \r
 722                 T minval_, maxval_;\r
 723                 cudaSafeCall(cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 724                 cudaSafeCall(cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 725                 *minval = minval_;\r
 726                 *maxval = maxval_;\r
 727 \r
 728                 uint minloc_, maxloc_;\r
 729                 cudaSafeCall(cudaMemcpy(&minloc_, minloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 730                 cudaSafeCall(cudaMemcpy(&maxloc_, maxloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 731                 minloc[1] = minloc_ / src.cols; minloc[0] = minloc_ - minloc[1] * src.cols;\r
 732                 maxloc[1] = maxloc_ / src.cols; maxloc[0] = maxloc_ - maxloc[1] * src.cols;\r
 733             }\r
 734 \r
 735             template void minMaxLocMaskMultipassCaller<uchar>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 736             template void minMaxLocMaskMultipassCaller<char>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 737             template void minMaxLocMaskMultipassCaller<ushort>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 738             template void minMaxLocMaskMultipassCaller<short>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 739             template void minMaxLocMaskMultipassCaller<int>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 740             template void minMaxLocMaskMultipassCaller<float>(const DevMem2Db, const PtrStepb, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 741 \r
 742 \r
 743             template <typename T>\r
 744             void minMaxLocMultipassCaller(const DevMem2Db src, double* minval, double* maxval, \r
 745                                           int minloc[2], int maxloc[2], PtrStepb valbuf, PtrStepb locbuf)\r
 746             {\r
 747                 dim3 threads, grid;\r
 748                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 749                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 750 \r
 751                 T* minval_buf = (T*)valbuf.ptr(0);\r
 752                 T* maxval_buf = (T*)valbuf.ptr(1);\r
 753                 uint* minloc_buf = (uint*)locbuf.ptr(0);\r
 754                 uint* maxloc_buf = (uint*)locbuf.ptr(1);\r
 755 \r
 756                 minMaxLocKernel<256, T, MaskTrue><<<grid, threads>>>(src, MaskTrue(), minval_buf, maxval_buf, \r
 757                                                                      minloc_buf, maxloc_buf);\r
 758                 cudaSafeCall( cudaGetLastError() );\r
 759                 minMaxLocPass2Kernel<256, T><<<1, 256>>>(minval_buf, maxval_buf, minloc_buf, maxloc_buf, grid.x * grid.y);\r
 760                 cudaSafeCall( cudaGetLastError() );\r
 761 \r
 762                 cudaSafeCall( cudaDeviceSynchronize() );\r
 763 \r
 764                 T minval_, maxval_;\r
 765                 cudaSafeCall(cudaMemcpy(&minval_, minval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 766                 cudaSafeCall(cudaMemcpy(&maxval_, maxval_buf, sizeof(T), cudaMemcpyDeviceToHost));\r
 767                 *minval = minval_;\r
 768                 *maxval = maxval_;\r
 769 \r
 770                 uint minloc_, maxloc_;\r
 771                 cudaSafeCall(cudaMemcpy(&minloc_, minloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 772                 cudaSafeCall(cudaMemcpy(&maxloc_, maxloc_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 773                 minloc[1] = minloc_ / src.cols; minloc[0] = minloc_ - minloc[1] * src.cols;\r
 774                 maxloc[1] = maxloc_ / src.cols; maxloc[0] = maxloc_ - maxloc[1] * src.cols;\r
 775             }\r
 776 \r
 777             template void minMaxLocMultipassCaller<uchar>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 778             template void minMaxLocMultipassCaller<char>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 779             template void minMaxLocMultipassCaller<ushort>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 780             template void minMaxLocMultipassCaller<short>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 781             template void minMaxLocMultipassCaller<int>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 782             template void minMaxLocMultipassCaller<float>(const DevMem2Db, double*, double*, int[2], int[2], PtrStepb, PtrStepb);\r
 783         } // namespace minmaxloc\r
 784 \r
 785         //////////////////////////////////////////////////////////////////////////////////////////////////////////\r
 786         // countNonZero\r
 787 \r
 788         namespace countnonzero \r
 789         {\r
 790             __constant__ int ctwidth;\r
 791             __constant__ int ctheight;\r
 792 \r
 793             __device__ uint blocks_finished = 0;\r
 794 \r
 795             void estimateThreadCfg(int cols, int rows, dim3& threads, dim3& grid)\r
 796             {\r
 797                 threads = dim3(32, 8);\r
 798                 grid = dim3(divUp(cols, threads.x * 8), divUp(rows, threads.y * 32));\r
 799                 grid.x = std::min(grid.x, threads.x);\r
 800                 grid.y = std::min(grid.y, threads.y);\r
 801             }\r
 802 \r
 803 \r
 804             void getBufSizeRequired(int cols, int rows, int& bufcols, int& bufrows)\r
 805             {\r
 806                 dim3 threads, grid;\r
 807                 estimateThreadCfg(cols, rows, threads, grid);\r
 808                 bufcols = grid.x * grid.y * sizeof(int);\r
 809                 bufrows = 1;\r
 810             }\r
 811 \r
 812 \r
 813             void setKernelConsts(int cols, int rows, const dim3& threads, const dim3& grid)\r
 814             {        \r
 815                 int twidth = divUp(divUp(cols, grid.x), threads.x);\r
 816                 int theight = divUp(divUp(rows, grid.y), threads.y);\r
 817                 cudaSafeCall(cudaMemcpyToSymbol(ctwidth, &twidth, sizeof(twidth))); \r
 818                 cudaSafeCall(cudaMemcpyToSymbol(ctheight, &theight, sizeof(theight))); \r
 819             }\r
 820 \r
 821 \r
 822             template <int nthreads, typename T>\r
 823             __global__ void countNonZeroKernel(const DevMem2Db src, volatile uint* count)\r
 824             {\r
 825                 __shared__ uint scount[nthreads];\r
 826 \r
 827                 uint x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
 828                 uint y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
 829                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 830 \r
 831                         uint cnt = 0;\r
 832                 for (uint y = 0; y < ctheight && y0 + y * blockDim.y < src.rows; ++y)\r
 833                 {\r
 834                     const T* ptr = (const T*)src.ptr(y0 + y * blockDim.y);\r
 835                     for (uint x = 0; x < ctwidth && x0 + x * blockDim.x < src.cols; ++x)\r
 836                                         cnt += ptr[x0 + x * blockDim.x] != 0;\r
 837                         }\r
 838 \r
 839                         scount[tid] = cnt;\r
 840                         __syncthreads();\r
 841 \r
 842                 sumInSmem<nthreads, uint>(scount, tid);\r
 843 \r
 844             #if __CUDA_ARCH__ >= 110\r
 845                         __shared__ bool is_last;\r
 846 \r
 847                         if (tid == 0)\r
 848                         {\r
 849                                 count[blockIdx.y * gridDim.x + blockIdx.x] = scount[0];\r
 850                                 __threadfence();\r
 851 \r
 852                                 uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
 853                                 is_last = ticket == gridDim.x * gridDim.y - 1;\r
 854                         }\r
 855 \r
 856                         __syncthreads();\r
 857 \r
 858                         if (is_last)\r
 859                         {\r
 860                     scount[tid] = tid < gridDim.x * gridDim.y ? count[tid] : 0;\r
 861                     __syncthreads();\r
 862 \r
 863                                 sumInSmem<nthreads, uint>(scount, tid);\r
 864 \r
 865                                 if (tid == 0) \r
 866                     {\r
 867                         count[0] = scount[0];\r
 868                         blocks_finished = 0;\r
 869                     }\r
 870                         }\r
 871             #else\r
 872                 if (tid == 0) count[blockIdx.y * gridDim.x + blockIdx.x] = scount[0];\r
 873             #endif\r
 874             }\r
 875 \r
 876            \r
 877             template <typename T>\r
 878             int countNonZeroCaller(const DevMem2Db src, PtrStepb buf)\r
 879             {\r
 880                 dim3 threads, grid;\r
 881                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 882                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 883 \r
 884                 uint* count_buf = (uint*)buf.ptr(0);\r
 885 \r
 886                 countNonZeroKernel<256, T><<<grid, threads>>>(src, count_buf);\r
 887                 cudaSafeCall( cudaGetLastError() );\r
 888 \r
 889                 cudaSafeCall( cudaDeviceSynchronize() );\r
 890 \r
 891                 uint count;\r
 892                 cudaSafeCall(cudaMemcpy(&count, count_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 893                 \r
 894                 return count;\r
 895             }  \r
 896 \r
 897             template int countNonZeroCaller<uchar>(const DevMem2Db, PtrStepb);\r
 898             template int countNonZeroCaller<char>(const DevMem2Db, PtrStepb);\r
 899             template int countNonZeroCaller<ushort>(const DevMem2Db, PtrStepb);\r
 900             template int countNonZeroCaller<short>(const DevMem2Db, PtrStepb);\r
 901             template int countNonZeroCaller<int>(const DevMem2Db, PtrStepb);\r
 902             template int countNonZeroCaller<float>(const DevMem2Db, PtrStepb);\r
 903             template int countNonZeroCaller<double>(const DevMem2Db, PtrStepb);\r
 904 \r
 905 \r
 906             template <int nthreads, typename T>\r
 907             __global__ void countNonZeroPass2Kernel(uint* count, int size)\r
 908             {\r
 909                 __shared__ uint scount[nthreads];\r
 910                 uint tid = threadIdx.y * blockDim.x + threadIdx.x;\r
 911 \r
 912                 scount[tid] = tid < size ? count[tid] : 0;\r
 913                 __syncthreads();\r
 914 \r
 915                 sumInSmem<nthreads, uint>(scount, tid);\r
 916 \r
 917                 if (tid == 0) \r
 918                     count[0] = scount[0];\r
 919             }\r
 920 \r
 921 \r
 922             template <typename T>\r
 923             int countNonZeroMultipassCaller(const DevMem2Db src, PtrStepb buf)\r
 924             {\r
 925                 dim3 threads, grid;\r
 926                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
 927                 setKernelConsts(src.cols, src.rows, threads, grid);\r
 928 \r
 929                 uint* count_buf = (uint*)buf.ptr(0);\r
 930 \r
 931                 countNonZeroKernel<256, T><<<grid, threads>>>(src, count_buf);\r
 932                 cudaSafeCall( cudaGetLastError() );\r
 933                 countNonZeroPass2Kernel<256, T><<<1, 256>>>(count_buf, grid.x * grid.y);\r
 934                 cudaSafeCall( cudaGetLastError() );\r
 935 \r
 936                 cudaSafeCall( cudaDeviceSynchronize() );\r
 937 \r
 938                 uint count;\r
 939                 cudaSafeCall(cudaMemcpy(&count, count_buf, sizeof(int), cudaMemcpyDeviceToHost));\r
 940                 \r
 941                 return count;\r
 942             }  \r
 943 \r
 944             template int countNonZeroMultipassCaller<uchar>(const DevMem2Db, PtrStepb);\r
 945             template int countNonZeroMultipassCaller<char>(const DevMem2Db, PtrStepb);\r
 946             template int countNonZeroMultipassCaller<ushort>(const DevMem2Db, PtrStepb);\r
 947             template int countNonZeroMultipassCaller<short>(const DevMem2Db, PtrStepb);\r
 948             template int countNonZeroMultipassCaller<int>(const DevMem2Db, PtrStepb);\r
 949             template int countNonZeroMultipassCaller<float>(const DevMem2Db, PtrStepb);\r
 950 \r
 951         } // namespace countnonzero\r
 952 \r
 953 \r
 954         //////////////////////////////////////////////////////////////////////////\r
 955         // Sum\r
 956 \r
 957         namespace sum\r
 958         {\r
 959             template <typename T> struct SumType {};\r
 960             template <> struct SumType<uchar> { typedef uint R; };\r
 961             template <> struct SumType<char> { typedef int R; };\r
 962             template <> struct SumType<ushort> { typedef uint R; };\r
 963             template <> struct SumType<short> { typedef int R; };\r
 964             template <> struct SumType<int> { typedef int R; };\r
 965             template <> struct SumType<float> { typedef float R; };\r
 966             template <> struct SumType<double> { typedef double R; };\r
 967 \r
 968             template <typename R> \r
 969             struct IdentityOp { static __device__ __forceinline__ R call(R x) { return x; } };\r
 970 \r
 971             template <typename R> \r
 972             struct AbsOp { static __device__ __forceinline__ R call(R x) { return ::abs(x); } };\r
 973 \r
 974             template <>\r
 975             struct AbsOp<uint> { static __device__ __forceinline__ uint call(uint x) { return x; } };\r
 976 \r
 977             template <typename R> \r
 978             struct SqrOp { static __device__ __forceinline__ R call(R x) { return x * x; } };\r
 979 \r
 980             __constant__ int ctwidth;\r
 981             __constant__ int ctheight;\r
 982             __device__ uint blocks_finished = 0;\r
 983 \r
 984             const int threads_x = 32;\r
 985             const int threads_y = 8;\r
 986 \r
 987             void estimateThreadCfg(int cols, int rows, dim3& threads, dim3& grid)\r
 988             {\r
 989                 threads = dim3(threads_x, threads_y);\r
 990                 grid = dim3(divUp(cols, threads.x * threads.y), \r
 991                             divUp(rows, threads.y * threads.x));\r
 992                 grid.x = std::min(grid.x, threads.x);\r
 993                 grid.y = std::min(grid.y, threads.y);\r
 994             }\r
 995 \r
 996 \r
 997             void getBufSizeRequired(int cols, int rows, int cn, int& bufcols, int& bufrows)\r
 998             {\r
 999                 dim3 threads, grid;\r
1000                 estimateThreadCfg(cols, rows, threads, grid);\r
1001                 bufcols = grid.x * grid.y * sizeof(double) * cn;\r
1002                 bufrows = 1;\r
1003             }\r
1004 \r
1005 \r
1006             void setKernelConsts(int cols, int rows, const dim3& threads, const dim3& grid)\r
1007             {        \r
1008                 int twidth = divUp(divUp(cols, grid.x), threads.x);\r
1009                 int theight = divUp(divUp(rows, grid.y), threads.y);\r
1010                 cudaSafeCall(cudaMemcpyToSymbol(ctwidth, &twidth, sizeof(twidth))); \r
1011                 cudaSafeCall(cudaMemcpyToSymbol(ctheight, &theight, sizeof(theight))); \r
1012             }\r
1013 \r
1014             template <typename T, typename R, typename Op, int nthreads>\r
1015             __global__ void sumKernel(const DevMem2Db src, R* result)\r
1016             {\r
1017                 __shared__ R smem[nthreads];\r
1018 \r
1019                 const int x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
1020                 const int y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
1021                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1022                 const int bid = blockIdx.y * gridDim.x + blockIdx.x;\r
1023 \r
1024                 R sum = 0;\r
1025                 for (int y = 0; y < ctheight && y0 + y * blockDim.y < src.rows; ++y)\r
1026                 {\r
1027                     const T* ptr = (const T*)src.ptr(y0 + y * blockDim.y);\r
1028                     for (int x = 0; x < ctwidth && x0 + x * blockDim.x < src.cols; ++x)\r
1029                         sum += Op::call(ptr[x0 + x * blockDim.x]);\r
1030                 }\r
1031 \r
1032                 smem[tid] = sum;\r
1033                 __syncthreads();\r
1034 \r
1035                 sumInSmem<nthreads, R>(smem, tid);\r
1036 \r
1037             #if __CUDA_ARCH__ >= 110\r
1038                 __shared__ bool is_last;\r
1039 \r
1040                 if (tid == 0)\r
1041                 {\r
1042                     result[bid] = smem[0];\r
1043                     __threadfence();\r
1044 \r
1045                     uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
1046                     is_last = (ticket == gridDim.x * gridDim.y - 1);\r
1047                 }\r
1048 \r
1049                 __syncthreads();\r
1050 \r
1051                 if (is_last)\r
1052                 {\r
1053                     smem[tid] = tid < gridDim.x * gridDim.y ? result[tid] : 0;\r
1054                     __syncthreads();\r
1055 \r
1056                     sumInSmem<nthreads, R>(smem, tid);\r
1057 \r
1058                     if (tid == 0) \r
1059                     {\r
1060                         result[0] = smem[0];\r
1061                         blocks_finished = 0;\r
1062                     }\r
1063                 }\r
1064             #else\r
1065                 if (tid == 0) result[bid] = smem[0];\r
1066             #endif\r
1067             }\r
1068 \r
1069 \r
1070             template <typename T, typename R, int nthreads>\r
1071             __global__ void sumPass2Kernel(R* result, int size)\r
1072             {\r
1073                 __shared__ R smem[nthreads];\r
1074                 int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1075 \r
1076                 smem[tid] = tid < size ? result[tid] : 0;\r
1077                 __syncthreads();\r
1078 \r
1079                 sumInSmem<nthreads, R>(smem, tid);\r
1080 \r
1081                 if (tid == 0) \r
1082                     result[0] = smem[0];\r
1083             }\r
1084 \r
1085 \r
1086             template <typename T, typename R, typename Op, int nthreads>\r
1087             __global__ void sumKernel_C2(const DevMem2Db src, typename TypeVec<R, 2>::vec_type* result)\r
1088             {\r
1089                 typedef typename TypeVec<T, 2>::vec_type SrcType;\r
1090                 typedef typename TypeVec<R, 2>::vec_type DstType;\r
1091 \r
1092                 __shared__ R smem[nthreads * 2];\r
1093 \r
1094                 const int x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
1095                 const int y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
1096                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1097                 const int bid = blockIdx.y * gridDim.x + blockIdx.x;\r
1098 \r
1099                 SrcType val;\r
1100                 DstType sum = VecTraits<DstType>::all(0);\r
1101                 for (int y = 0; y < ctheight && y0 + y * blockDim.y < src.rows; ++y)\r
1102                 {\r
1103                     const SrcType* ptr = (const SrcType*)src.ptr(y0 + y * blockDim.y);\r
1104                     for (int x = 0; x < ctwidth && x0 + x * blockDim.x < src.cols; ++x)\r
1105                     {\r
1106                         val = ptr[x0 + x * blockDim.x];\r
1107                         sum = sum + VecTraits<DstType>::make(Op::call(val.x), Op::call(val.y));\r
1108                     }\r
1109                 }\r
1110 \r
1111                 smem[tid] = sum.x;\r
1112                 smem[tid + nthreads] = sum.y;\r
1113                 __syncthreads();\r
1114 \r
1115                 sumInSmem<nthreads, R>(smem, tid);\r
1116                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1117 \r
1118             #if __CUDA_ARCH__ >= 110\r
1119                 __shared__ bool is_last;\r
1120 \r
1121                 if (tid == 0)\r
1122                 {\r
1123                     DstType res;\r
1124                     res.x = smem[0];\r
1125                     res.y = smem[nthreads];\r
1126                     result[bid] = res;\r
1127                     __threadfence();\r
1128 \r
1129                     uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
1130                     is_last = (ticket == gridDim.x * gridDim.y - 1);\r
1131                 }\r
1132 \r
1133                 __syncthreads();\r
1134 \r
1135                 if (is_last)\r
1136                 {\r
1137                     DstType res = tid < gridDim.x * gridDim.y ? result[tid] : VecTraits<DstType>::all(0);\r
1138                     smem[tid] = res.x;\r
1139                     smem[tid + nthreads] = res.y;\r
1140                     __syncthreads();\r
1141 \r
1142                     sumInSmem<nthreads, R>(smem, tid);\r
1143                     sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1144 \r
1145                     if (tid == 0) \r
1146                     {\r
1147                         res.x = smem[0];\r
1148                         res.y = smem[nthreads];\r
1149                         result[0] = res;\r
1150                         blocks_finished = 0;\r
1151                     }\r
1152                 }\r
1153             #else\r
1154                 if (tid == 0) \r
1155                 {\r
1156                     DstType res;\r
1157                     res.x = smem[0];\r
1158                     res.y = smem[nthreads];\r
1159                     result[bid] = res;\r
1160                 }\r
1161             #endif\r
1162             }\r
1163 \r
1164 \r
1165             template <typename T, typename R, int nthreads>\r
1166             __global__ void sumPass2Kernel_C2(typename TypeVec<R, 2>::vec_type* result, int size)\r
1167             {\r
1168                 typedef typename TypeVec<R, 2>::vec_type DstType;\r
1169 \r
1170                 __shared__ R smem[nthreads * 2];\r
1171 \r
1172                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1173 \r
1174                 DstType res = tid < size ? result[tid] : VecTraits<DstType>::all(0);\r
1175                 smem[tid] = res.x;\r
1176                 smem[tid + nthreads] = res.y;\r
1177                 __syncthreads();\r
1178 \r
1179                 sumInSmem<nthreads, R>(smem, tid);\r
1180                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1181 \r
1182                 if (tid == 0) \r
1183                 {\r
1184                     res.x = smem[0];\r
1185                     res.y = smem[nthreads];\r
1186                     result[0] = res;\r
1187                 }\r
1188             }\r
1189 \r
1190 \r
1191             template <typename T, typename R, typename Op, int nthreads>\r
1192             __global__ void sumKernel_C3(const DevMem2Db src, typename TypeVec<R, 3>::vec_type* result)\r
1193             {\r
1194                 typedef typename TypeVec<T, 3>::vec_type SrcType;\r
1195                 typedef typename TypeVec<R, 3>::vec_type DstType;\r
1196 \r
1197                 __shared__ R smem[nthreads * 3];\r
1198 \r
1199                 const int x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
1200                 const int y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
1201                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1202                 const int bid = blockIdx.y * gridDim.x + blockIdx.x;\r
1203 \r
1204                 SrcType val;\r
1205                 DstType sum = VecTraits<DstType>::all(0);\r
1206                 for (int y = 0; y < ctheight && y0 + y * blockDim.y < src.rows; ++y)\r
1207                 {\r
1208                     const SrcType* ptr = (const SrcType*)src.ptr(y0 + y * blockDim.y);\r
1209                     for (int x = 0; x < ctwidth && x0 + x * blockDim.x < src.cols; ++x)\r
1210                     {\r
1211                         val = ptr[x0 + x * blockDim.x];\r
1212                         sum = sum + VecTraits<DstType>::make(Op::call(val.x), Op::call(val.y), Op::call(val.z));\r
1213                     }\r
1214                 }\r
1215 \r
1216                 smem[tid] = sum.x;\r
1217                 smem[tid + nthreads] = sum.y;\r
1218                 smem[tid + 2 * nthreads] = sum.z;\r
1219                 __syncthreads();\r
1220 \r
1221                 sumInSmem<nthreads, R>(smem, tid);\r
1222                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1223                 sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1224 \r
1225             #if __CUDA_ARCH__ >= 110\r
1226                 __shared__ bool is_last;\r
1227 \r
1228                 if (tid == 0)\r
1229                 {\r
1230                     DstType res;\r
1231                     res.x = smem[0];\r
1232                     res.y = smem[nthreads];\r
1233                     res.z = smem[2 * nthreads];\r
1234                     result[bid] = res;\r
1235                     __threadfence();\r
1236 \r
1237                     uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
1238                     is_last = (ticket == gridDim.x * gridDim.y - 1);\r
1239                 }\r
1240 \r
1241                 __syncthreads();\r
1242 \r
1243                 if (is_last)\r
1244                 {\r
1245                     DstType res = tid < gridDim.x * gridDim.y ? result[tid] : VecTraits<DstType>::all(0);\r
1246                     smem[tid] = res.x;\r
1247                     smem[tid + nthreads] = res.y;\r
1248                     smem[tid + 2 * nthreads] = res.z;\r
1249                     __syncthreads();\r
1250 \r
1251                     sumInSmem<nthreads, R>(smem, tid);\r
1252                     sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1253                     sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1254 \r
1255                     if (tid == 0) \r
1256                     {\r
1257                         res.x = smem[0];\r
1258                         res.y = smem[nthreads];\r
1259                         res.z = smem[2 * nthreads];\r
1260                         result[0] = res;\r
1261                         blocks_finished = 0;\r
1262                     }\r
1263                 }\r
1264             #else\r
1265                 if (tid == 0) \r
1266                 {\r
1267                     DstType res;\r
1268                     res.x = smem[0];\r
1269                     res.y = smem[nthreads];\r
1270                     res.z = smem[2 * nthreads];\r
1271                     result[bid] = res;\r
1272                 }\r
1273             #endif\r
1274             }\r
1275 \r
1276 \r
1277             template <typename T, typename R, int nthreads>\r
1278             __global__ void sumPass2Kernel_C3(typename TypeVec<R, 3>::vec_type* result, int size)\r
1279             {\r
1280                 typedef typename TypeVec<R, 3>::vec_type DstType;\r
1281 \r
1282                 __shared__ R smem[nthreads * 3];\r
1283 \r
1284                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1285 \r
1286                 DstType res = tid < size ? result[tid] : VecTraits<DstType>::all(0);\r
1287                 smem[tid] = res.x;\r
1288                 smem[tid + nthreads] = res.y;\r
1289                 smem[tid + 2 * nthreads] = res.z;\r
1290                 __syncthreads();\r
1291 \r
1292                 sumInSmem<nthreads, R>(smem, tid);\r
1293                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1294                 sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1295 \r
1296                 if (tid == 0) \r
1297                 {\r
1298                     res.x = smem[0];\r
1299                     res.y = smem[nthreads];\r
1300                     res.z = smem[2 * nthreads];\r
1301                     result[0] = res;\r
1302                 }\r
1303             }\r
1304 \r
1305             template <typename T, typename R, typename Op, int nthreads>\r
1306             __global__ void sumKernel_C4(const DevMem2Db src, typename TypeVec<R, 4>::vec_type* result)\r
1307             {\r
1308                 typedef typename TypeVec<T, 4>::vec_type SrcType;\r
1309                 typedef typename TypeVec<R, 4>::vec_type DstType;\r
1310 \r
1311                 __shared__ R smem[nthreads * 4];\r
1312 \r
1313                 const int x0 = blockIdx.x * blockDim.x * ctwidth + threadIdx.x;\r
1314                 const int y0 = blockIdx.y * blockDim.y * ctheight + threadIdx.y;\r
1315                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1316                 const int bid = blockIdx.y * gridDim.x + blockIdx.x;\r
1317 \r
1318                 SrcType val;\r
1319                 DstType sum = VecTraits<DstType>::all(0);\r
1320                 for (int y = 0; y < ctheight && y0 + y * blockDim.y < src.rows; ++y)\r
1321                 {\r
1322                     const SrcType* ptr = (const SrcType*)src.ptr(y0 + y * blockDim.y);\r
1323                     for (int x = 0; x < ctwidth && x0 + x * blockDim.x < src.cols; ++x)\r
1324                     {\r
1325                         val = ptr[x0 + x * blockDim.x];\r
1326                         sum = sum + VecTraits<DstType>::make(Op::call(val.x), Op::call(val.y), \r
1327                                                              Op::call(val.z), Op::call(val.w));\r
1328                     }\r
1329                 }\r
1330 \r
1331                 smem[tid] = sum.x;\r
1332                 smem[tid + nthreads] = sum.y;\r
1333                 smem[tid + 2 * nthreads] = sum.z;\r
1334                 smem[tid + 3 * nthreads] = sum.w;\r
1335                 __syncthreads();\r
1336 \r
1337                 sumInSmem<nthreads, R>(smem, tid);\r
1338                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1339                 sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1340                 sumInSmem<nthreads, R>(smem + 3 * nthreads, tid);\r
1341 \r
1342             #if __CUDA_ARCH__ >= 110\r
1343                 __shared__ bool is_last;\r
1344 \r
1345                 if (tid == 0)\r
1346                 {\r
1347                     DstType res;\r
1348                     res.x = smem[0];\r
1349                     res.y = smem[nthreads];\r
1350                     res.z = smem[2 * nthreads];\r
1351                     res.w = smem[3 * nthreads];\r
1352                     result[bid] = res;\r
1353                     __threadfence();\r
1354 \r
1355                     uint ticket = atomicInc(&blocks_finished, gridDim.x * gridDim.y);\r
1356                     is_last = (ticket == gridDim.x * gridDim.y - 1);\r
1357                 }\r
1358 \r
1359                 __syncthreads();\r
1360 \r
1361                 if (is_last)\r
1362                 {\r
1363                     DstType res = tid < gridDim.x * gridDim.y ? result[tid] : VecTraits<DstType>::all(0);\r
1364                     smem[tid] = res.x;\r
1365                     smem[tid + nthreads] = res.y;\r
1366                     smem[tid + 2 * nthreads] = res.z;\r
1367                     smem[tid + 3 * nthreads] = res.w;\r
1368                     __syncthreads();\r
1369 \r
1370                     sumInSmem<nthreads, R>(smem, tid);\r
1371                     sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1372                     sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1373                     sumInSmem<nthreads, R>(smem + 3 * nthreads, tid);\r
1374 \r
1375                     if (tid == 0) \r
1376                     {\r
1377                         res.x = smem[0];\r
1378                         res.y = smem[nthreads];\r
1379                         res.z = smem[2 * nthreads];\r
1380                         res.w = smem[3 * nthreads];\r
1381                         result[0] = res;\r
1382                         blocks_finished = 0;\r
1383                     }\r
1384                 }\r
1385             #else\r
1386                 if (tid == 0) \r
1387                 {\r
1388                     DstType res;\r
1389                     res.x = smem[0];\r
1390                     res.y = smem[nthreads];\r
1391                     res.z = smem[2 * nthreads];\r
1392                     res.w = smem[3 * nthreads];\r
1393                     result[bid] = res;\r
1394                 }\r
1395             #endif\r
1396             }\r
1397 \r
1398 \r
1399             template <typename T, typename R, int nthreads>\r
1400             __global__ void sumPass2Kernel_C4(typename TypeVec<R, 4>::vec_type* result, int size)\r
1401             {\r
1402                 typedef typename TypeVec<R, 4>::vec_type DstType;\r
1403 \r
1404                 __shared__ R smem[nthreads * 4];\r
1405 \r
1406                 const int tid = threadIdx.y * blockDim.x + threadIdx.x;\r
1407 \r
1408                 DstType res = tid < size ? result[tid] : VecTraits<DstType>::all(0);\r
1409                 smem[tid] = res.x;\r
1410                 smem[tid + nthreads] = res.y;\r
1411                 smem[tid + 2 * nthreads] = res.z;\r
1412                 smem[tid + 3 * nthreads] = res.w;\r
1413                 __syncthreads();\r
1414 \r
1415                 sumInSmem<nthreads, R>(smem, tid);\r
1416                 sumInSmem<nthreads, R>(smem + nthreads, tid);\r
1417                 sumInSmem<nthreads, R>(smem + 2 * nthreads, tid);\r
1418                 sumInSmem<nthreads, R>(smem + 3 * nthreads, tid);\r
1419 \r
1420                 if (tid == 0) \r
1421                 {\r
1422                     res.x = smem[0];\r
1423                     res.y = smem[nthreads];\r
1424                     res.z = smem[2 * nthreads];\r
1425                     res.w = smem[3 * nthreads];\r
1426                     result[0] = res;\r
1427                 }\r
1428             }\r
1429 \r
1430             template <typename T>\r
1431             void sumMultipassCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1432             {\r
1433                 typedef typename SumType<T>::R R;\r
1434 \r
1435                 dim3 threads, grid;\r
1436                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1437                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1438 \r
1439                 switch (cn)\r
1440                 {\r
1441                 case 1:\r
1442                     sumKernel<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1443                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1444                     cudaSafeCall( cudaGetLastError() );\r
1445 \r
1446                     sumPass2Kernel<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1447                             (typename TypeVec<R, 1>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1448                     cudaSafeCall( cudaGetLastError() );\r
1449 \r
1450                     break;\r
1451                 case 2:\r
1452                     sumKernel_C2<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1453                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1454                     cudaSafeCall( cudaGetLastError() );\r
1455 \r
1456                     sumPass2Kernel_C2<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1457                             (typename TypeVec<R, 2>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1458                     cudaSafeCall( cudaGetLastError() );\r
1459 \r
1460                     break;\r
1461                 case 3:\r
1462                     sumKernel_C3<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1463                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1464                     cudaSafeCall( cudaGetLastError() );\r
1465 \r
1466                     sumPass2Kernel_C3<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1467                             (typename TypeVec<R, 3>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1468                     cudaSafeCall( cudaGetLastError() );\r
1469 \r
1470                     break;\r
1471                 case 4:\r
1472                     sumKernel_C4<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1473                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1474                     cudaSafeCall( cudaGetLastError() );\r
1475 \r
1476                     sumPass2Kernel_C4<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1477                             (typename TypeVec<R, 4>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1478                     cudaSafeCall( cudaGetLastError() );\r
1479 \r
1480                     break;\r
1481                 }\r
1482                 cudaSafeCall( cudaDeviceSynchronize() );\r
1483 \r
1484                 R result[4] = {0, 0, 0, 0};\r
1485                 cudaSafeCall(cudaMemcpy(&result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1486 \r
1487                 sum[0] = result[0];\r
1488                 sum[1] = result[1];\r
1489                 sum[2] = result[2];\r
1490                 sum[3] = result[3];\r
1491             }  \r
1492 \r
1493             template void sumMultipassCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1494             template void sumMultipassCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1495             template void sumMultipassCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1496             template void sumMultipassCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1497             template void sumMultipassCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1498             template void sumMultipassCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1499 \r
1500 \r
1501             template <typename T>\r
1502             void sumCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1503             {\r
1504                 typedef typename SumType<T>::R R;\r
1505 \r
1506                 dim3 threads, grid;\r
1507                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1508                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1509 \r
1510                 switch (cn)\r
1511                 {\r
1512                 case 1:\r
1513                     sumKernel<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1514                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1515                     break;\r
1516                 case 2:\r
1517                     sumKernel_C2<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1518                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1519                     break;\r
1520                 case 3:\r
1521                     sumKernel_C3<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1522                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1523                     break;\r
1524                 case 4:\r
1525                     sumKernel_C4<T, R, IdentityOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1526                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1527                     break;\r
1528                 }\r
1529                 cudaSafeCall( cudaGetLastError() );\r
1530 \r
1531                 cudaSafeCall( cudaDeviceSynchronize() );\r
1532 \r
1533                 R result[4] = {0, 0, 0, 0};\r
1534                 cudaSafeCall(cudaMemcpy(&result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1535 \r
1536                 sum[0] = result[0];\r
1537                 sum[1] = result[1];\r
1538                 sum[2] = result[2];\r
1539                 sum[3] = result[3];\r
1540             }  \r
1541 \r
1542             template void sumCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1543             template void sumCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1544             template void sumCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1545             template void sumCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1546             template void sumCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1547             template void sumCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1548 \r
1549 \r
1550             template <typename T>\r
1551             void absSumMultipassCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1552             {\r
1553                 typedef typename SumType<T>::R R;\r
1554 \r
1555                 dim3 threads, grid;\r
1556                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1557                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1558 \r
1559                 switch (cn)\r
1560                 {\r
1561                 case 1:\r
1562                     sumKernel<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1563                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1564                     cudaSafeCall( cudaGetLastError() );\r
1565 \r
1566                     sumPass2Kernel<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1567                             (typename TypeVec<R, 1>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1568                     cudaSafeCall( cudaGetLastError() );\r
1569 \r
1570                     break;\r
1571                 case 2:\r
1572                     sumKernel_C2<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1573                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1574                     cudaSafeCall( cudaGetLastError() );\r
1575 \r
1576                     sumPass2Kernel_C2<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1577                             (typename TypeVec<R, 2>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1578                     cudaSafeCall( cudaGetLastError() );\r
1579 \r
1580                     break;\r
1581                 case 3:\r
1582                     sumKernel_C3<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1583                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1584                     cudaSafeCall( cudaGetLastError() );\r
1585 \r
1586                     sumPass2Kernel_C3<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1587                             (typename TypeVec<R, 3>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1588                     cudaSafeCall( cudaGetLastError() );\r
1589 \r
1590                     break;\r
1591                 case 4:\r
1592                     sumKernel_C4<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1593                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1594                     cudaSafeCall( cudaGetLastError() );\r
1595 \r
1596                     sumPass2Kernel_C4<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1597                             (typename TypeVec<R, 4>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1598                     cudaSafeCall( cudaGetLastError() );\r
1599 \r
1600                     break;\r
1601                 }\r
1602                 cudaSafeCall( cudaDeviceSynchronize() );\r
1603 \r
1604                 R result[4] = {0, 0, 0, 0};\r
1605                 cudaSafeCall(cudaMemcpy(result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1606 \r
1607                 sum[0] = result[0];\r
1608                 sum[1] = result[1];\r
1609                 sum[2] = result[2];\r
1610                 sum[3] = result[3];\r
1611             }  \r
1612 \r
1613             template void absSumMultipassCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1614             template void absSumMultipassCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1615             template void absSumMultipassCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1616             template void absSumMultipassCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1617             template void absSumMultipassCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1618             template void absSumMultipassCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1619 \r
1620 \r
1621             template <typename T>\r
1622             void absSumCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1623             {\r
1624                 typedef typename SumType<T>::R R;\r
1625 \r
1626                 dim3 threads, grid;\r
1627                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1628                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1629 \r
1630                 switch (cn)\r
1631                 {\r
1632                 case 1:\r
1633                     sumKernel<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1634                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1635                     break;\r
1636                 case 2:\r
1637                     sumKernel_C2<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1638                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1639                     break;\r
1640                 case 3:\r
1641                     sumKernel_C3<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1642                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1643                     break;\r
1644                 case 4:\r
1645                     sumKernel_C4<T, R, AbsOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1646                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1647                     break;\r
1648                 }\r
1649                 cudaSafeCall( cudaGetLastError() );\r
1650 \r
1651                 cudaSafeCall( cudaDeviceSynchronize() );\r
1652 \r
1653                 R result[4] = {0, 0, 0, 0};\r
1654                 cudaSafeCall(cudaMemcpy(result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1655 \r
1656                 sum[0] = result[0];\r
1657                 sum[1] = result[1];\r
1658                 sum[2] = result[2];\r
1659                 sum[3] = result[3];\r
1660             }\r
1661 \r
1662             template void absSumCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1663             template void absSumCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1664             template void absSumCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1665             template void absSumCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1666             template void absSumCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1667             template void absSumCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1668 \r
1669 \r
1670             template <typename T>\r
1671             void sqrSumMultipassCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1672             {\r
1673                 typedef typename SumType<T>::R R;\r
1674 \r
1675                 dim3 threads, grid;\r
1676                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1677                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1678 \r
1679                 switch (cn)\r
1680                 {\r
1681                 case 1:\r
1682                     sumKernel<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1683                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1684                     cudaSafeCall( cudaGetLastError() );\r
1685 \r
1686                     sumPass2Kernel<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1687                             (typename TypeVec<R, 1>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1688                     cudaSafeCall( cudaGetLastError() );\r
1689 \r
1690                     break;\r
1691                 case 2:\r
1692                     sumKernel_C2<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1693                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1694                     cudaSafeCall( cudaGetLastError() );\r
1695 \r
1696                     sumPass2Kernel_C2<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1697                             (typename TypeVec<R, 2>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1698                     cudaSafeCall( cudaGetLastError() );\r
1699 \r
1700                     break;\r
1701                 case 3:\r
1702                     sumKernel_C3<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1703                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1704                     cudaSafeCall( cudaGetLastError() );\r
1705 \r
1706                     sumPass2Kernel_C3<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1707                             (typename TypeVec<R, 3>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1708                     cudaSafeCall( cudaGetLastError() );\r
1709 \r
1710                     break;\r
1711                 case 4:\r
1712                     sumKernel_C4<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1713                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1714                     cudaSafeCall( cudaGetLastError() );\r
1715 \r
1716                     sumPass2Kernel_C4<T, R, threads_x * threads_y><<<1, threads_x * threads_y>>>(\r
1717                             (typename TypeVec<R, 4>::vec_type*)buf.ptr(0), grid.x * grid.y);\r
1718                     cudaSafeCall( cudaGetLastError() );\r
1719 \r
1720                     break;\r
1721                 }\r
1722                 cudaSafeCall( cudaDeviceSynchronize() );\r
1723 \r
1724                 R result[4] = {0, 0, 0, 0};\r
1725                 cudaSafeCall(cudaMemcpy(result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1726 \r
1727                 sum[0] = result[0];\r
1728                 sum[1] = result[1];\r
1729                 sum[2] = result[2];\r
1730                 sum[3] = result[3];\r
1731             }  \r
1732 \r
1733             template void sqrSumMultipassCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1734             template void sqrSumMultipassCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1735             template void sqrSumMultipassCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1736             template void sqrSumMultipassCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1737             template void sqrSumMultipassCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1738             template void sqrSumMultipassCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1739 \r
1740 \r
1741             template <typename T>\r
1742             void sqrSumCaller(const DevMem2Db src, PtrStepb buf, double* sum, int cn)\r
1743             {\r
1744                 typedef double R;\r
1745 \r
1746                 dim3 threads, grid;\r
1747                 estimateThreadCfg(src.cols, src.rows, threads, grid);\r
1748                 setKernelConsts(src.cols, src.rows, threads, grid);\r
1749 \r
1750                 switch (cn)\r
1751                 {\r
1752                 case 1:\r
1753                     sumKernel<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1754                             src, (typename TypeVec<R, 1>::vec_type*)buf.ptr(0));\r
1755                     break;\r
1756                 case 2:\r
1757                     sumKernel_C2<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1758                             src, (typename TypeVec<R, 2>::vec_type*)buf.ptr(0));\r
1759                     break;\r
1760                 case 3:\r
1761                     sumKernel_C3<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1762                             src, (typename TypeVec<R, 3>::vec_type*)buf.ptr(0));\r
1763                     break;\r
1764                 case 4:\r
1765                     sumKernel_C4<T, R, SqrOp<R>, threads_x * threads_y><<<grid, threads>>>(\r
1766                             src, (typename TypeVec<R, 4>::vec_type*)buf.ptr(0));\r
1767                     break;\r
1768                 }\r
1769                 cudaSafeCall( cudaGetLastError() );\r
1770 \r
1771                 cudaSafeCall( cudaDeviceSynchronize() );\r
1772 \r
1773                 R result[4] = {0, 0, 0, 0};\r
1774                 cudaSafeCall(cudaMemcpy(result, buf.ptr(0), sizeof(R) * cn, cudaMemcpyDeviceToHost));\r
1775 \r
1776                 sum[0] = result[0];\r
1777                 sum[1] = result[1];\r
1778                 sum[2] = result[2];\r
1779                 sum[3] = result[3];\r
1780             }\r
1781 \r
1782             template void sqrSumCaller<uchar>(const DevMem2Db, PtrStepb, double*, int);\r
1783             template void sqrSumCaller<char>(const DevMem2Db, PtrStepb, double*, int);\r
1784             template void sqrSumCaller<ushort>(const DevMem2Db, PtrStepb, double*, int);\r
1785             template void sqrSumCaller<short>(const DevMem2Db, PtrStepb, double*, int);\r
1786             template void sqrSumCaller<int>(const DevMem2Db, PtrStepb, double*, int);\r
1787             template void sqrSumCaller<float>(const DevMem2Db, PtrStepb, double*, int);\r
1788         } // namespace sum\r
1789 \r
1790         //////////////////////////////////////////////////////////////////////////////\r
1791         // reduce\r
1792 \r
1793         template <typename S> struct SumReductor\r
1794         {\r
1795             __device__ __forceinline__ S startValue() const\r
1796             {\r
1797                 return 0;\r
1798             }\r
1799 \r
1800             __device__ __forceinline__ SumReductor(const SumReductor& other){}\r
1801             __device__ __forceinline__ SumReductor(){}\r
1802 \r
1803             __device__ __forceinline__ S operator ()(volatile S a, volatile S b) const\r
1804             {\r
1805                 return a + b;\r
1806             }\r
1807 \r
1808             __device__ __forceinline__ S result(S r, double) const\r
1809             {\r
1810                 return r;\r
1811             }\r
1812         };\r
1813 \r
1814         template <typename S> struct AvgReductor\r
1815         {\r
1816             __device__ __forceinline__ S startValue() const\r
1817             {\r
1818                 return 0;\r
1819             }\r
1820 \r
1821             __device__ __forceinline__ AvgReductor(const AvgReductor& other){}\r
1822             __device__ __forceinline__ AvgReductor(){}\r
1823 \r
1824             __device__ __forceinline__ S operator ()(volatile S a, volatile S b) const\r
1825             {\r
1826                 return a + b;\r
1827             }\r
1828 \r
1829             __device__ __forceinline__ double result(S r, double sz) const\r
1830             {\r
1831                 return r / sz;\r
1832             }\r
1833         };\r
1834 \r
1835         template <typename S> struct MinReductor\r
1836         {\r
1837             __device__ __forceinline__ S startValue() const\r
1838             {\r
1839                 return numeric_limits<S>::max();\r
1840             }\r
1841 \r
1842             __device__ __forceinline__ MinReductor(const MinReductor& other){}\r
1843             __device__ __forceinline__ MinReductor(){}\r
1844 \r
1845             template <typename T> __device__ __forceinline__ T operator ()(volatile T a, volatile T b) const\r
1846             {\r
1847                 return saturate_cast<T>(::min(a, b));\r
1848             }\r
1849             __device__ __forceinline__ float operator ()(volatile float a, volatile float b) const\r
1850             {\r
1851                 return ::fmin(a, b);\r
1852             }\r
1853 \r
1854             __device__ __forceinline__ S result(S r, double) const\r
1855             {\r
1856                 return r;\r
1857             }\r
1858         };\r
1859 \r
1860         template <typename S> struct MaxReductor\r
1861         {\r
1862             __device__ __forceinline__ S startValue() const\r
1863             {\r
1864                 return numeric_limits<S>::min();\r
1865             }\r
1866 \r
1867             __device__ __forceinline__ MaxReductor(const MaxReductor& other){}\r
1868             __device__ __forceinline__ MaxReductor(){}\r
1869 \r
1870             template <typename T> __device__ __forceinline__ int operator ()(volatile T a, volatile T b) const\r
1871             {\r
1872                 return ::max(a, b);\r
1873             }\r
1874             __device__ __forceinline__ float operator ()(volatile float a, volatile float b) const\r
1875             {\r
1876                 return ::fmax(a, b);\r
1877             }\r
1878 \r
1879             __device__ __forceinline__ S result(S r, double) const\r
1880             {\r
1881                 return r;\r
1882             }\r
1883         };\r
1884 \r
1885         template <class Op, typename T, typename S, typename D> __global__ void reduceRows(const DevMem2D_<T> src, D* dst, const Op op)\r
1886         {\r
1887             __shared__ S smem[16 * 16];\r
1888 \r
1889             const int x = blockIdx.x * 16 + threadIdx.x;\r
1890 \r
1891             S myVal = op.startValue();\r
1892 \r
1893             if (x < src.cols)\r
1894             {\r
1895                 for (int y = threadIdx.y; y < src.rows; y += 16)\r
1896                     myVal = op(myVal, src.ptr(y)[x]);\r
1897             }        \r
1898 \r
1899             smem[threadIdx.x * 16 + threadIdx.y] = myVal;\r
1900             __syncthreads();\r
1901 \r
1902             if (threadIdx.x < 8)\r
1903             {\r
1904                 volatile S* srow = smem + threadIdx.y * 16;\r
1905                 srow[threadIdx.x] = op(srow[threadIdx.x], srow[threadIdx.x + 8]);\r
1906                 srow[threadIdx.x] = op(srow[threadIdx.x], srow[threadIdx.x + 4]);\r
1907                 srow[threadIdx.x] = op(srow[threadIdx.x], srow[threadIdx.x + 2]);\r
1908                 srow[threadIdx.x] = op(srow[threadIdx.x], srow[threadIdx.x + 1]);\r
1909             }\r
1910             __syncthreads();\r
1911 \r
1912             if (threadIdx.y == 0 && x < src.cols)\r
1913                 dst[x] = saturate_cast<D>(op.result(smem[threadIdx.x * 16], src.rows));\r
1914         }\r
1915 \r
1916         template <template <typename> class Op, typename T, typename S, typename D> void reduceRows_caller(const DevMem2D_<T>& src, DevMem2D_<D> dst, cudaStream_t stream)\r
1917         {\r
1918             const dim3 block(16, 16);\r
1919             const dim3 grid(divUp(src.cols, block.x));\r
1920 \r
1921             Op<S> op;\r
1922             reduceRows<Op<S>, T, S, D><<<grid, block, 0, stream>>>(src, dst.data, op);\r
1923             cudaSafeCall( cudaGetLastError() );\r
1924 \r
1925             if (stream == 0)\r
1926                 cudaSafeCall( cudaDeviceSynchronize() );\r
1927 \r
1928         }\r
1929 \r
1930         template <typename T, typename S, typename D> void reduceRows_gpu(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream)\r
1931         {\r
1932             typedef void (*caller_t)(const DevMem2D_<T>& src, DevMem2D_<D> dst, cudaStream_t stream);\r
1933 \r
1934             static const caller_t callers[] = \r
1935             {\r
1936                 reduceRows_caller<SumReductor, T, S, D>, \r
1937                 reduceRows_caller<AvgReductor, T, S, D>, \r
1938                 reduceRows_caller<MaxReductor, T, S, D>, \r
1939                 reduceRows_caller<MinReductor, T, S, D>\r
1940             };\r
1941 \r
1942             callers[reduceOp](static_cast< DevMem2D_<T> >(src), static_cast< DevMem2D_<D> >(dst), stream);\r
1943         }\r
1944 \r
1945         template void reduceRows_gpu<uchar, int, uchar>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1946         template void reduceRows_gpu<uchar, int, int>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1947         template void reduceRows_gpu<uchar, int, float>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);  \r
1948 \r
1949         template void reduceRows_gpu<ushort, int, ushort>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1950         template void reduceRows_gpu<ushort, int, int>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1951         template void reduceRows_gpu<ushort, int, float>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream); \r
1952 \r
1953         template void reduceRows_gpu<short, int, short>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1954         template void reduceRows_gpu<short, int, int>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1955         template void reduceRows_gpu<short, int, float>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream); \r
1956 \r
1957         template void reduceRows_gpu<int, int, int>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1958         template void reduceRows_gpu<int, int, float>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1959 \r
1960         template void reduceRows_gpu<float, float, float>(const DevMem2Db& src, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
1961 \r
1962 \r
1963 \r
1964         template <int cn, class Op, typename T, typename S, typename D> __global__ void reduceCols(const DevMem2D_<T> src, D* dst, const Op op)\r
1965         {\r
1966             __shared__ S smem[256 * cn];\r
1967 \r
1968             const int y = blockIdx.x;\r
1969 \r
1970             const T* src_row = src.ptr(y);\r
1971 \r
1972             S myVal[cn];\r
1973 \r
1974             #pragma unroll\r
1975             for (int c = 0; c < cn; ++c)\r
1976                 myVal[c] = op.startValue();\r
1977 \r
1978         #if __CUDA_ARCH__ >= 200\r
1979 \r
1980             // For cc >= 2.0 prefer L1 cache\r
1981             for (int x = threadIdx.x; x < src.cols; x += 256)\r
1982             {\r
1983                 #pragma unroll\r
1984                 for (int c = 0; c < cn; ++c)\r
1985                     myVal[c] = op(myVal[c], src_row[x * cn + c]);\r
1986             }\r
1987 \r
1988         #else // __CUDA_ARCH__ >= 200\r
1989 \r
1990             // For older arch use shared memory for cache\r
1991             for (int x = 0; x < src.cols; x += 256)\r
1992             {\r
1993                 #pragma unroll\r
1994                 for (int c = 0; c < cn; ++c)\r
1995                 {\r
1996                     smem[c * 256 + threadIdx.x] = op.startValue();\r
1997                     const int load_x = x * cn + c * 256 + threadIdx.x;\r
1998                     if (load_x < src.cols * cn)\r
1999                         smem[c * 256 + threadIdx.x] = src_row[load_x];\r
2000                 }\r
2001                 __syncthreads();\r
2002 \r
2003                 #pragma unroll\r
2004                 for (int c = 0; c < cn; ++c)\r
2005                     myVal[c] = op(myVal[c], smem[threadIdx.x * cn + c]);\r
2006                 __syncthreads();\r
2007             }\r
2008 \r
2009         #endif // __CUDA_ARCH__ >= 200\r
2010 \r
2011             #pragma unroll\r
2012             for (int c = 0; c < cn; ++c)\r
2013                 smem[c * 256 + threadIdx.x] = myVal[c];\r
2014             __syncthreads();\r
2015 \r
2016             if (threadIdx.x < 128)\r
2017             {\r
2018                 #pragma unroll\r
2019                 for (int c = 0; c < cn; ++c)\r
2020                     smem[c * 256 + threadIdx.x] = op(smem[c * 256 + threadIdx.x], smem[c * 256 + threadIdx.x + 128]);\r
2021             }\r
2022             __syncthreads();\r
2023 \r
2024             if (threadIdx.x < 64)\r
2025             {\r
2026                 #pragma unroll\r
2027                 for (int c = 0; c < cn; ++c)\r
2028                     smem[c * 256 + threadIdx.x] = op(smem[c * 256 + threadIdx.x], smem[c * 256 + threadIdx.x + 64]);\r
2029             }\r
2030             __syncthreads();\r
2031 \r
2032             volatile S* sdata = smem;\r
2033 \r
2034             if (threadIdx.x < 32)\r
2035             {\r
2036                 #pragma unroll\r
2037                 for (int c = 0; c < cn; ++c)\r
2038                 {\r
2039                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 32]);\r
2040                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 16]);\r
2041                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 8]);\r
2042                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 4]);\r
2043                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 2]);\r
2044                     sdata[c * 256 + threadIdx.x] = op(sdata[c * 256 + threadIdx.x], sdata[c * 256 + threadIdx.x + 1]);\r
2045                 }\r
2046             }\r
2047             __syncthreads();\r
2048 \r
2049             if (threadIdx.x < cn)\r
2050                 dst[y * cn + threadIdx.x] = saturate_cast<D>(op.result(smem[threadIdx.x * 256], src.cols));\r
2051         }\r
2052 \r
2053         template <int cn, template <typename> class Op, typename T, typename S, typename D> void reduceCols_caller(const DevMem2D_<T>& src, DevMem2D_<D> dst, cudaStream_t stream)\r
2054         {\r
2055             const dim3 block(256);\r
2056             const dim3 grid(src.rows);\r
2057 \r
2058             Op<S> op;\r
2059             reduceCols<cn, Op<S>, T, S, D><<<grid, block, 0, stream>>>(src, dst.data, op);\r
2060             cudaSafeCall( cudaGetLastError() );\r
2061 \r
2062             if (stream == 0)\r
2063                 cudaSafeCall( cudaDeviceSynchronize() );\r
2064 \r
2065         }\r
2066 \r
2067         template <typename T, typename S, typename D> void reduceCols_gpu(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream)\r
2068         {\r
2069             typedef void (*caller_t)(const DevMem2D_<T>& src, DevMem2D_<D> dst, cudaStream_t stream);\r
2070 \r
2071             static const caller_t callers[4][4] = \r
2072             {\r
2073                 {reduceCols_caller<1, SumReductor, T, S, D>, reduceCols_caller<1, AvgReductor, T, S, D>, reduceCols_caller<1, MaxReductor, T, S, D>, reduceCols_caller<1, MinReductor, T, S, D>},\r
2074                 {reduceCols_caller<2, SumReductor, T, S, D>, reduceCols_caller<2, AvgReductor, T, S, D>, reduceCols_caller<2, MaxReductor, T, S, D>, reduceCols_caller<2, MinReductor, T, S, D>},\r
2075                 {reduceCols_caller<3, SumReductor, T, S, D>, reduceCols_caller<3, AvgReductor, T, S, D>, reduceCols_caller<3, MaxReductor, T, S, D>, reduceCols_caller<3, MinReductor, T, S, D>},\r
2076                 {reduceCols_caller<4, SumReductor, T, S, D>, reduceCols_caller<4, AvgReductor, T, S, D>, reduceCols_caller<4, MaxReductor, T, S, D>, reduceCols_caller<4, MinReductor, T, S, D>},\r
2077             };\r
2078 \r
2079             callers[cn - 1][reduceOp](static_cast< DevMem2D_<T> >(src), static_cast< DevMem2D_<D> >(dst), stream);\r
2080         }\r
2081 \r
2082         template void reduceCols_gpu<uchar, int, uchar>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2083         template void reduceCols_gpu<uchar, int, int>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2084         template void reduceCols_gpu<uchar, int, float>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2085 \r
2086         template void reduceCols_gpu<ushort, int, ushort>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream); \r
2087         template void reduceCols_gpu<ushort, int, int>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);                  \r
2088         template void reduceCols_gpu<ushort, int, float>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2089 \r
2090         template void reduceCols_gpu<short, int, short>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);  \r
2091         template void reduceCols_gpu<short, int, int>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);                  \r
2092         template void reduceCols_gpu<short, int, float>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);  \r
2093 \r
2094         template void reduceCols_gpu<int, int, int>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);                  \r
2095         template void reduceCols_gpu<int, int, float>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2096 \r
2097         template void reduceCols_gpu<float, float, float>(const DevMem2Db& src, int cn, const DevMem2Db& dst, int reduceOp, cudaStream_t stream);\r
2098     } // namespace mattrix_reductions\r
2099 }}} // namespace cv { namespace gpu { namespace device\r