modules/gpu/src/cuda/bf_knnmatch.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 bpied warranties, including, but not limited to, the bpied\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/vec_distance.hpp"\r
  46 \r
  47 using namespace cv::gpu;\r
  48 using namespace cv::gpu::device;\r
  49 \r
  50 namespace cv { namespace gpu { namespace bf_knnmatch\r
  51 {\r
  52     ///////////////////////////////////////////////////////////////////////////////\r
  53     // Reduction\r
  54 \r
  55     template <int BLOCK_SIZE> \r
  56     __device__ void findBestMatch(float& bestDistance1, float& bestDistance2, \r
  57                                   int& bestTrainIdx1, int& bestTrainIdx2, \r
  58                                   float* s_distance, int* s_trainIdx)\r
  59     {\r
  60         float myBestDistance1 = numeric_limits<float>::max(); \r
  61         float myBestDistance2 = numeric_limits<float>::max();\r
  62         int myBestTrainIdx1 = -1;\r
  63         int myBestTrainIdx2 = -1;\r
  64 \r
  65         s_distance += threadIdx.y * BLOCK_SIZE;\r
  66         s_trainIdx += threadIdx.y * BLOCK_SIZE;\r
  67 \r
  68         s_distance[threadIdx.x] = bestDistance1;\r
  69         s_trainIdx[threadIdx.x] = bestTrainIdx1;\r
  70 \r
  71         __syncthreads();\r
  72 \r
  73         if (threadIdx.x == 0)\r
  74         {\r
  75             #pragma unroll\r
  76             for (int i = 0; i < BLOCK_SIZE; ++i)\r
  77             {\r
  78                 float val = s_distance[i];\r
  79 \r
  80                 if (val < myBestDistance1)\r
  81                 {\r
  82                     myBestDistance2 = myBestDistance1;\r
  83                     myBestTrainIdx2 = myBestTrainIdx1;\r
  84 \r
  85                     myBestDistance1 = val;\r
  86                     myBestTrainIdx1 = s_trainIdx[i];\r
  87                 }\r
  88                 else if (val < myBestDistance2)\r
  89                 {\r
  90                     myBestDistance2 = val;\r
  91                     myBestTrainIdx2 = s_trainIdx[i];\r
  92                 }\r
  93             }\r
  94         }\r
  95 \r
  96         __syncthreads();\r
  97 \r
  98         s_distance[threadIdx.x] = bestDistance2;\r
  99         s_trainIdx[threadIdx.x] = bestTrainIdx2;\r
 100 \r
 101         __syncthreads();\r
 102 \r
 103         if (threadIdx.x == 0)\r
 104         {\r
 105             #pragma unroll\r
 106             for (int i = 0; i < BLOCK_SIZE; ++i)\r
 107             {\r
 108                 float val = s_distance[i];\r
 109 \r
 110                 if (val < myBestDistance2)\r
 111                 {\r
 112                     myBestDistance2 = val;\r
 113                     myBestTrainIdx2 = s_trainIdx[i];\r
 114                 }\r
 115             }\r
 116         }\r
 117 \r
 118         bestDistance1 = myBestDistance1;\r
 119         bestDistance2 = myBestDistance2;\r
 120 \r
 121         bestTrainIdx1 = myBestTrainIdx1;\r
 122         bestTrainIdx2 = myBestTrainIdx2;\r
 123     }\r
 124 \r
 125     template <int BLOCK_SIZE> \r
 126     __device__ void findBestMatch(float& bestDistance1, float& bestDistance2, \r
 127                                    int& bestTrainIdx1, int& bestTrainIdx2, \r
 128                                    int& bestImgIdx1, int& bestImgIdx2, \r
 129                                    float* s_distance, int* s_trainIdx, int* s_imgIdx)\r
 130     {\r
 131         float myBestDistance1 = numeric_limits<float>::max(); \r
 132         float myBestDistance2 = numeric_limits<float>::max();\r
 133         int myBestTrainIdx1 = -1;\r
 134         int myBestTrainIdx2 = -1;\r
 135         int myBestImgIdx1 = -1;\r
 136         int myBestImgIdx2 = -1;\r
 137 \r
 138         s_distance += threadIdx.y * BLOCK_SIZE;\r
 139         s_trainIdx += threadIdx.y * BLOCK_SIZE;\r
 140         s_imgIdx   += threadIdx.y * BLOCK_SIZE;\r
 141 \r
 142         s_distance[threadIdx.x] = bestDistance1;\r
 143         s_trainIdx[threadIdx.x] = bestTrainIdx1;\r
 144         s_imgIdx[threadIdx.x]   = bestImgIdx1;\r
 145 \r
 146         __syncthreads();\r
 147 \r
 148         if (threadIdx.x == 0)\r
 149         {\r
 150             #pragma unroll\r
 151             for (int i = 0; i < BLOCK_SIZE; ++i)\r
 152             {\r
 153                 float val = s_distance[i];\r
 154 \r
 155                 if (val < myBestDistance1)\r
 156                 {\r
 157                     myBestDistance2 = myBestDistance1;\r
 158                     myBestTrainIdx2 = myBestTrainIdx1;\r
 159                     myBestImgIdx2   = myBestImgIdx1;\r
 160 \r
 161                     myBestDistance1 = val;\r
 162                     myBestTrainIdx1 = s_trainIdx[i];\r
 163                     myBestImgIdx1   = s_imgIdx[i];\r
 164                 }\r
 165                 else if (val < myBestDistance2)\r
 166                 {\r
 167                     myBestDistance2 = val;\r
 168                     myBestTrainIdx2 = s_trainIdx[i];\r
 169                     myBestImgIdx2   = s_imgIdx[i];\r
 170                 }\r
 171             }\r
 172         }\r
 173 \r
 174         __syncthreads();\r
 175 \r
 176         s_distance[threadIdx.x] = bestDistance2;\r
 177         s_trainIdx[threadIdx.x] = bestTrainIdx2;\r
 178         s_imgIdx[threadIdx.x]   = bestImgIdx2;\r
 179 \r
 180         __syncthreads();\r
 181 \r
 182         if (threadIdx.x == 0)\r
 183         {\r
 184             #pragma unroll\r
 185             for (int i = 0; i < BLOCK_SIZE; ++i)\r
 186             {\r
 187                 float val = s_distance[i];\r
 188 \r
 189                 if (val < myBestDistance2)\r
 190                 {\r
 191                     myBestDistance2 = val;\r
 192                     myBestTrainIdx2 = s_trainIdx[i];\r
 193                     myBestImgIdx2   = s_imgIdx[i];\r
 194                 }\r
 195             }\r
 196         }\r
 197 \r
 198         bestDistance1 = myBestDistance1;\r
 199         bestDistance2 = myBestDistance2;\r
 200 \r
 201         bestTrainIdx1 = myBestTrainIdx1;\r
 202         bestTrainIdx2 = myBestTrainIdx2;\r
 203 \r
 204         bestImgIdx1 = myBestImgIdx1;\r
 205         bestImgIdx2 = myBestImgIdx2;\r
 206     }\r
 207 \r
 208     ///////////////////////////////////////////////////////////////////////////////\r
 209     // Match Unrolled Cached\r
 210 \r
 211     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename T, typename U> \r
 212     __device__ void loadQueryToSmem(int queryIdx, const DevMem2D_<T>& query, U* s_query)\r
 213     {\r
 214         #pragma unroll\r
 215         for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; ++i)\r
 216         {\r
 217             const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 218             s_query[threadIdx.y * MAX_DESC_LEN + loadX] = loadX < query.cols ? query.ptr(min(queryIdx, query.rows - 1))[loadX] : 0;\r
 219         }\r
 220     }\r
 221 \r
 222     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 223     __device__ void loopUnrolledCached(int queryIdx, const DevMem2D_<T>& query, int imgIdx, const DevMem2D_<T>& train, const Mask& mask, \r
 224                                        typename Dist::value_type* s_query, typename Dist::value_type* s_train, \r
 225                                        float& bestDistance1, float& bestDistance2, \r
 226                                        int& bestTrainIdx1, int& bestTrainIdx2, \r
 227                                        int& bestImgIdx1, int& bestImgIdx2)\r
 228     {\r
 229         for (int t = 0, endt = (train.rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; ++t)\r
 230         {\r
 231             Dist dist;\r
 232 \r
 233             #pragma unroll\r
 234             for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; ++i)\r
 235             {\r
 236                 const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 237 \r
 238                 s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = loadX < train.cols ? train.ptr(min(t * BLOCK_SIZE + threadIdx.y, train.rows - 1))[loadX] : 0;\r
 239 \r
 240                 __syncthreads();\r
 241 \r
 242                 #pragma unroll\r
 243                 for (int j = 0; j < BLOCK_SIZE; ++j)\r
 244                     dist.reduceIter(s_query[threadIdx.y * MAX_DESC_LEN + i * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + threadIdx.x]);\r
 245 \r
 246                 __syncthreads();\r
 247             }\r
 248 \r
 249             typename Dist::result_type distVal = dist;\r
 250 \r
 251             const int trainIdx = t * BLOCK_SIZE + threadIdx.x;\r
 252 \r
 253             if (queryIdx < query.rows && trainIdx < train.rows && mask(queryIdx, trainIdx))\r
 254             {\r
 255                 if (distVal < bestDistance1)\r
 256                 {\r
 257                     bestImgIdx2   = bestImgIdx1;\r
 258                     bestDistance2 = bestDistance1;\r
 259                     bestTrainIdx2 = bestTrainIdx1;\r
 260 \r
 261                     bestImgIdx1   = imgIdx;\r
 262                     bestDistance1 = distVal;\r
 263                     bestTrainIdx1 = trainIdx;\r
 264                 }\r
 265                 else if (distVal < bestDistance2)\r
 266                 {\r
 267                     bestImgIdx2   = imgIdx;\r
 268                     bestDistance2 = distVal;\r
 269                     bestTrainIdx2 = trainIdx;\r
 270                 }\r
 271             }\r
 272         }\r
 273     }\r
 274 \r
 275     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 276     __global__ void matchUnrolledCached(const DevMem2D_<T> query, const DevMem2D_<T> train, const Mask mask, int2* bestTrainIdx, float2* bestDistance)\r
 277     {\r
 278         extern __shared__ int smem[];\r
 279 \r
 280         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 281 \r
 282         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 283         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * MAX_DESC_LEN);\r
 284 \r
 285         loadQueryToSmem<BLOCK_SIZE, MAX_DESC_LEN>(queryIdx, query, s_query);\r
 286 \r
 287         float myBestDistance1 = numeric_limits<float>::max();\r
 288         float myBestDistance2 = numeric_limits<float>::max();\r
 289         int myBestTrainIdx1 = -1;\r
 290         int myBestTrainIdx2 = -1;\r
 291 \r
 292         loopUnrolledCached<BLOCK_SIZE, MAX_DESC_LEN, Dist>(queryIdx, query, 0, train, mask, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestTrainIdx1, myBestTrainIdx2);\r
 293 \r
 294         __syncthreads();\r
 295 \r
 296         float* s_distance = (float*)(smem);\r
 297         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 298 \r
 299         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, s_distance, s_trainIdx);\r
 300 \r
 301         if (queryIdx < query.rows && threadIdx.x == 0)\r
 302         {\r
 303             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 304             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 305         }\r
 306     }\r
 307 \r
 308     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 309     void matchUnrolledCached(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, \r
 310                              const DevMem2D_<int2>& trainIdx, const DevMem2D_<float2>& distance, \r
 311                              cudaStream_t stream)\r
 312     {\r
 313         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 314         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 315 \r
 316         const size_t smemSize = (BLOCK_SIZE * (MAX_DESC_LEN >= BLOCK_SIZE ? MAX_DESC_LEN : BLOCK_SIZE) + BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 317 \r
 318         matchUnrolledCached<BLOCK_SIZE, MAX_DESC_LEN, Dist><<<grid, block, smemSize, stream>>>(query, train, mask, trainIdx.data, distance.data);\r
 319         cudaSafeCall( cudaGetLastError() );\r
 320 \r
 321         if (stream == 0)\r
 322             cudaSafeCall( cudaDeviceSynchronize() );\r
 323     }\r
 324 \r
 325     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 326     __global__ void matchUnrolledCached(const DevMem2D_<T> query, const DevMem2D_<T>* trains, int n, const Mask mask, int2* bestTrainIdx, int2* bestImgIdx, float2* bestDistance)\r
 327     {\r
 328         extern __shared__ int smem[];\r
 329 \r
 330         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 331 \r
 332         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 333         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * MAX_DESC_LEN);\r
 334 \r
 335         loadQueryToSmem<BLOCK_SIZE, MAX_DESC_LEN>(queryIdx, query, s_query);\r
 336 \r
 337         float myBestDistance1 = numeric_limits<float>::max();\r
 338         float myBestDistance2 = numeric_limits<float>::max();\r
 339         int myBestTrainIdx1 = -1;\r
 340         int myBestTrainIdx2 = -1;\r
 341         int myBestImgIdx1 = -1;\r
 342         int myBestImgIdx2 = -1;\r
 343 \r
 344         Mask m = mask;\r
 345 \r
 346         for (int imgIdx = 0; imgIdx < n; ++imgIdx)\r
 347         {\r
 348             const DevMem2D_<T> train = trains[imgIdx];\r
 349             m.next();\r
 350             loopUnrolledCached<BLOCK_SIZE, MAX_DESC_LEN, Dist>(queryIdx, query, imgIdx, train, m, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2);\r
 351         }\r
 352 \r
 353         __syncthreads();\r
 354 \r
 355         float* s_distance = (float*)(smem);\r
 356         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 357         int* s_imgIdx = (int*)(smem + 2 * BLOCK_SIZE * BLOCK_SIZE);\r
 358 \r
 359         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2, s_distance, s_trainIdx, s_imgIdx);\r
 360 \r
 361         if (queryIdx < query.rows && threadIdx.x == 0)\r
 362         {\r
 363             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 364             bestImgIdx[queryIdx] = make_int2(myBestImgIdx1, myBestImgIdx2);\r
 365             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 366         }\r
 367     }\r
 368 \r
 369     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 370     void matchUnrolledCached(const DevMem2D_<T>& query, const DevMem2D_<T>* trains, int n, const Mask& mask, \r
 371                              const DevMem2D_<int2>& trainIdx, const DevMem2D_<int2>& imgIdx, const DevMem2D_<float2>& distance, \r
 372                              cudaStream_t stream)\r
 373     {\r
 374         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 375         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 376 \r
 377         const size_t smemSize = (BLOCK_SIZE * (MAX_DESC_LEN >= 2 * BLOCK_SIZE ? MAX_DESC_LEN : 2 * BLOCK_SIZE) + BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 378 \r
 379         matchUnrolledCached<BLOCK_SIZE, MAX_DESC_LEN, Dist><<<grid, block, smemSize, stream>>>(query, trains, n, mask, trainIdx.data, imgIdx.data, distance.data);\r
 380         cudaSafeCall( cudaGetLastError() );\r
 381 \r
 382         if (stream == 0)\r
 383             cudaSafeCall( cudaDeviceSynchronize() );\r
 384     }\r
 385 \r
 386     ///////////////////////////////////////////////////////////////////////////////\r
 387     // Match Unrolled\r
 388 \r
 389     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 390     __device__ void loopUnrolled(int queryIdx, const DevMem2D_<T>& query, int imgIdx, const DevMem2D_<T>& train, const Mask& mask, \r
 391                                  typename Dist::value_type* s_query, typename Dist::value_type* s_train, \r
 392                                  float& bestDistance1, float& bestDistance2, \r
 393                                  int& bestTrainIdx1, int& bestTrainIdx2, \r
 394                                  int& bestImgIdx1, int& bestImgIdx2)\r
 395     {\r
 396         for (int t = 0, endt = (train.rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; ++t)\r
 397         {\r
 398             Dist dist;\r
 399 \r
 400             #pragma unroll\r
 401             for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; ++i)\r
 402             {\r
 403                 const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 404 \r
 405                 if (loadX < query.cols)\r
 406                 {\r
 407                     s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = query.ptr(min(queryIdx, query.rows - 1))[loadX];\r
 408                     s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = train.ptr(min(t * BLOCK_SIZE + threadIdx.y, train.rows - 1))[loadX];\r
 409                 }\r
 410                 else\r
 411                 {                \r
 412                     s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = 0;\r
 413                     s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = 0;\r
 414                 }\r
 415 \r
 416                 __syncthreads();\r
 417 \r
 418                 #pragma unroll\r
 419                 for (int j = 0; j < BLOCK_SIZE; ++j)\r
 420                     dist.reduceIter(s_query[threadIdx.y * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + threadIdx.x]);\r
 421 \r
 422                 __syncthreads();\r
 423             }\r
 424 \r
 425             typename Dist::result_type distVal = dist;\r
 426 \r
 427             const int trainIdx = t * BLOCK_SIZE + threadIdx.x;\r
 428 \r
 429             if (queryIdx < query.rows && trainIdx < train.rows && mask(queryIdx, trainIdx))\r
 430             {\r
 431                 if (distVal < bestDistance1)\r
 432                 {\r
 433                     bestImgIdx2   = bestImgIdx1;\r
 434                     bestDistance2 = bestDistance1;\r
 435                     bestTrainIdx2 = bestTrainIdx1;\r
 436 \r
 437                     bestImgIdx1   = imgIdx;\r
 438                     bestDistance1 = distVal;\r
 439                     bestTrainIdx1 = trainIdx;\r
 440                 }\r
 441                 else if (distVal < bestDistance2)\r
 442                 {\r
 443                     bestImgIdx2   = imgIdx;\r
 444                     bestDistance2 = distVal;\r
 445                     bestTrainIdx2 = trainIdx;\r
 446                 }\r
 447             }\r
 448         }\r
 449     }\r
 450 \r
 451     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 452     __global__ void matchUnrolled(const DevMem2D_<T> query, const DevMem2D_<T> train, const Mask mask, int2* bestTrainIdx, float2* bestDistance)\r
 453     {\r
 454         extern __shared__ int smem[];\r
 455 \r
 456         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 457 \r
 458         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 459         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 460 \r
 461         float myBestDistance1 = numeric_limits<float>::max();\r
 462         float myBestDistance2 = numeric_limits<float>::max();\r
 463         int myBestTrainIdx1 = -1;\r
 464         int myBestTrainIdx2 = -1;\r
 465 \r
 466         loopUnrolled<BLOCK_SIZE, MAX_DESC_LEN, Dist>(queryIdx, query, 0, train, mask, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestTrainIdx1, myBestTrainIdx2);\r
 467 \r
 468         __syncthreads();\r
 469 \r
 470         float* s_distance = (float*)(smem);\r
 471         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 472 \r
 473         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, s_distance, s_trainIdx);\r
 474 \r
 475         if (queryIdx < query.rows && threadIdx.x == 0)\r
 476         {\r
 477             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 478             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 479         }\r
 480     }\r
 481 \r
 482     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 483     void matchUnrolled(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, \r
 484                        const DevMem2D_<int2>& trainIdx, const DevMem2D_<float2>& distance, \r
 485                        cudaStream_t stream)\r
 486     {\r
 487         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 488         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 489 \r
 490         const size_t smemSize = (2 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 491 \r
 492         matchUnrolled<BLOCK_SIZE, MAX_DESC_LEN, Dist><<<grid, block, smemSize, stream>>>(query, train, mask, trainIdx.data, distance.data);\r
 493         cudaSafeCall( cudaGetLastError() );\r
 494 \r
 495         if (stream == 0)\r
 496             cudaSafeCall( cudaDeviceSynchronize() );\r
 497     }\r
 498 \r
 499     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 500     __global__ void matchUnrolled(const DevMem2D_<T> query, const DevMem2D_<T>* trains, int n, const Mask mask, int2* bestTrainIdx, int2* bestImgIdx, float2* bestDistance)\r
 501     {\r
 502         extern __shared__ int smem[];\r
 503 \r
 504         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 505 \r
 506         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 507         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 508 \r
 509         float myBestDistance1 = numeric_limits<float>::max();\r
 510         float myBestDistance2 = numeric_limits<float>::max();\r
 511         int myBestTrainIdx1 = -1;\r
 512         int myBestTrainIdx2 = -1;\r
 513         int myBestImgIdx1 = -1;\r
 514         int myBestImgIdx2 = -1;\r
 515 \r
 516         Mask m = mask;\r
 517 \r
 518         for (int imgIdx = 0; imgIdx < n; ++imgIdx)\r
 519         {\r
 520             const DevMem2D_<T> train = trains[imgIdx];\r
 521             m.next();\r
 522             loopUnrolled<BLOCK_SIZE, MAX_DESC_LEN, Dist>(queryIdx, query, imgIdx, train, m, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2);\r
 523         }\r
 524 \r
 525         __syncthreads();\r
 526 \r
 527         float* s_distance = (float*)(smem);\r
 528         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 529         int* s_imgIdx = (int*)(smem + 2 * BLOCK_SIZE * BLOCK_SIZE);\r
 530 \r
 531         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2, s_distance, s_trainIdx, s_imgIdx);\r
 532 \r
 533         if (queryIdx < query.rows && threadIdx.x == 0)\r
 534         {\r
 535             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 536             bestImgIdx[queryIdx] = make_int2(myBestImgIdx1, myBestImgIdx2);\r
 537             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 538         }\r
 539     }\r
 540 \r
 541     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 542     void matchUnrolled(const DevMem2D_<T>& query, const DevMem2D_<T>* trains, int n, const Mask& mask, \r
 543                        const DevMem2D_<int2>& trainIdx, const DevMem2D_<int2>& imgIdx, const DevMem2D_<float2>& distance, \r
 544                        cudaStream_t stream)\r
 545     {\r
 546         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 547         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 548 \r
 549         const size_t smemSize = (3 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 550 \r
 551         matchUnrolled<BLOCK_SIZE, MAX_DESC_LEN, Dist><<<grid, block, smemSize, stream>>>(query, trains, n, mask, trainIdx.data, imgIdx.data, distance.data);\r
 552         cudaSafeCall( cudaGetLastError() );\r
 553 \r
 554         if (stream == 0)\r
 555             cudaSafeCall( cudaDeviceSynchronize() );\r
 556     }\r
 557 \r
 558     ///////////////////////////////////////////////////////////////////////////////\r
 559     // Match\r
 560 \r
 561     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 562     __device__ void loop(int queryIdx, const DevMem2D_<T>& query, int imgIdx, const DevMem2D_<T>& train, const Mask& mask, \r
 563                          typename Dist::value_type* s_query, typename Dist::value_type* s_train, \r
 564                          float& bestDistance1, float& bestDistance2, \r
 565                          int& bestTrainIdx1, int& bestTrainIdx2, \r
 566                          int& bestImgIdx1, int& bestImgIdx2)\r
 567     {\r
 568         for (int t = 0, endt = (train.rows + BLOCK_SIZE - 1) / BLOCK_SIZE; t < endt; ++t)\r
 569         {\r
 570             Dist dist;\r
 571 \r
 572             for (int i = 0, endi = (query.cols + BLOCK_SIZE - 1) / BLOCK_SIZE; i < endi; ++i)\r
 573             {\r
 574                 const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 575 \r
 576                 if (loadX < query.cols)\r
 577                 {\r
 578                     s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = query.ptr(min(queryIdx, query.rows - 1))[loadX];\r
 579                     s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = train.ptr(min(t * BLOCK_SIZE + threadIdx.y, train.rows - 1))[loadX];\r
 580                 }\r
 581                 else\r
 582                 {                \r
 583                     s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = 0;\r
 584                     s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = 0;\r
 585                 }\r
 586 \r
 587                 __syncthreads();\r
 588 \r
 589                 #pragma unroll\r
 590                 for (int j = 0; j < BLOCK_SIZE; ++j)\r
 591                     dist.reduceIter(s_query[threadIdx.y * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + threadIdx.x]);\r
 592 \r
 593                 __syncthreads();\r
 594             }\r
 595 \r
 596             typename Dist::result_type distVal = dist;\r
 597 \r
 598             const int trainIdx = t * BLOCK_SIZE + threadIdx.x;\r
 599 \r
 600             if (queryIdx < query.rows && trainIdx < train.rows && mask(queryIdx, trainIdx))\r
 601             {\r
 602                 if (distVal < bestDistance1)\r
 603                 {\r
 604                     bestImgIdx2   = bestImgIdx1;\r
 605                     bestDistance2 = bestDistance1;\r
 606                     bestTrainIdx2 = bestTrainIdx1;\r
 607 \r
 608                     bestImgIdx1   = imgIdx;\r
 609                     bestDistance1 = distVal;\r
 610                     bestTrainIdx1 = trainIdx;\r
 611                 }\r
 612                 else if (distVal < bestDistance2)\r
 613                 {\r
 614                     bestImgIdx2   = imgIdx;\r
 615                     bestDistance2 = distVal;\r
 616                     bestTrainIdx2 = trainIdx;\r
 617                 }\r
 618             }\r
 619         }\r
 620     }\r
 621 \r
 622     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 623     __global__ void match(const DevMem2D_<T> query, const DevMem2D_<T> train, const Mask mask, int2* bestTrainIdx, float2* bestDistance)\r
 624     {\r
 625         extern __shared__ int smem[];\r
 626 \r
 627         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 628 \r
 629         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 630         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 631 \r
 632         float myBestDistance1 = numeric_limits<float>::max();\r
 633         float myBestDistance2 = numeric_limits<float>::max();\r
 634         int myBestTrainIdx1 = -1;\r
 635         int myBestTrainIdx2 = -1;\r
 636 \r
 637         loop<BLOCK_SIZE, Dist>(queryIdx, query, 0, train, mask, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestTrainIdx1, myBestTrainIdx2);\r
 638 \r
 639         __syncthreads();\r
 640 \r
 641         float* s_distance = (float*)(smem);\r
 642         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 643 \r
 644         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, s_distance, s_trainIdx);\r
 645 \r
 646         if (queryIdx < query.rows && threadIdx.x == 0)\r
 647         {\r
 648             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 649             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 650         }\r
 651     }\r
 652 \r
 653     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 654     void match(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, \r
 655                const DevMem2D_<int2>& trainIdx, const DevMem2D_<float2>& distance, \r
 656                cudaStream_t stream)\r
 657     {\r
 658         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 659         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 660 \r
 661         const size_t smemSize = (2 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 662 \r
 663         match<BLOCK_SIZE, Dist><<<grid, block, smemSize, stream>>>(query, train, mask, trainIdx.data, distance.data);\r
 664         cudaSafeCall( cudaGetLastError() );\r
 665 \r
 666         if (stream == 0)\r
 667             cudaSafeCall( cudaDeviceSynchronize() );\r
 668     }\r
 669 \r
 670     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 671     __global__ void match(const DevMem2D_<T> query, const DevMem2D_<T>* trains, int n, const Mask mask, int2* bestTrainIdx, int2* bestImgIdx, float2* bestDistance)\r
 672     {\r
 673         extern __shared__ int smem[];\r
 674 \r
 675         const int queryIdx = blockIdx.x * BLOCK_SIZE + threadIdx.y;\r
 676 \r
 677         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 678         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 679 \r
 680         float myBestDistance1 = numeric_limits<float>::max();\r
 681         float myBestDistance2 = numeric_limits<float>::max();\r
 682         int myBestTrainIdx1 = -1;\r
 683         int myBestTrainIdx2 = -1;\r
 684         int myBestImgIdx1 = -1;\r
 685         int myBestImgIdx2 = -1;\r
 686 \r
 687         Mask m = mask;\r
 688 \r
 689         for (int imgIdx = 0; imgIdx < n; ++imgIdx)\r
 690         {\r
 691             const DevMem2D_<T> train = trains[imgIdx];\r
 692             m.next();\r
 693             loop<BLOCK_SIZE, Dist>(queryIdx, query, imgIdx, train, m, s_query, s_train, myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2);\r
 694         }\r
 695 \r
 696         __syncthreads();\r
 697 \r
 698         float* s_distance = (float*)(smem);\r
 699         int* s_trainIdx = (int*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 700         int* s_imgIdx = (int*)(smem + 2 * BLOCK_SIZE * BLOCK_SIZE);\r
 701 \r
 702         findBestMatch<BLOCK_SIZE>(myBestDistance1, myBestDistance2, myBestTrainIdx1, myBestTrainIdx2, myBestImgIdx1, myBestImgIdx2, s_distance, s_trainIdx, s_imgIdx);\r
 703 \r
 704         if (queryIdx < query.rows && threadIdx.x == 0)\r
 705         {\r
 706             bestTrainIdx[queryIdx] = make_int2(myBestTrainIdx1, myBestTrainIdx2);\r
 707             bestImgIdx[queryIdx] = make_int2(myBestImgIdx1, myBestImgIdx2);\r
 708             bestDistance[queryIdx] = make_float2(myBestDistance1, myBestDistance2);\r
 709         }\r
 710     }\r
 711 \r
 712     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 713     void match(const DevMem2D_<T>& query, const DevMem2D_<T>* trains, int n, const Mask& mask, \r
 714                const DevMem2D_<int2>& trainIdx, const DevMem2D_<int2>& imgIdx, const DevMem2D_<float2>& distance, \r
 715                cudaStream_t stream)\r
 716     {\r
 717         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 718         const dim3 grid(divUp(query.rows, BLOCK_SIZE));\r
 719 \r
 720         const size_t smemSize = (3 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 721 \r
 722         match<BLOCK_SIZE, Dist><<<grid, block, smemSize, stream>>>(query, trains, n, mask, trainIdx.data, imgIdx.data, distance.data);\r
 723         cudaSafeCall( cudaGetLastError() );\r
 724 \r
 725         if (stream == 0)\r
 726             cudaSafeCall( cudaDeviceSynchronize() );\r
 727     }\r
 728 \r
 729     ///////////////////////////////////////////////////////////////////////////////\r
 730     // knnMatch 2 dispatcher\r
 731 \r
 732     template <typename Dist, typename T, typename Mask> \r
 733     void match2Dispatcher(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, \r
 734                           const DevMem2Db& trainIdx, const DevMem2Db& distance, \r
 735                           int cc, cudaStream_t stream)\r
 736     {\r
 737         if (query.cols <= 64)\r
 738         {\r
 739             matchUnrolledCached<16, 64, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 740         }\r
 741         else if (query.cols <= 128)\r
 742         {\r
 743             matchUnrolledCached<16, 128, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 744         }\r
 745         /*else if (query.cols <= 256)\r
 746         {\r
 747             matchUnrolled<16, 256, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 748         }\r
 749         else if (query.cols <= 512)\r
 750         {            \r
 751             matchUnrolled<16, 512, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 752         }\r
 753         else if (query.cols <= 1024)\r
 754         {            \r
 755             matchUnrolled<16, 1024, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 756         }*/\r
 757         else\r
 758         {\r
 759             match<16, Dist>(query, train, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 760         }\r
 761     }\r
 762 \r
 763     template <typename Dist, typename T, typename Mask> \r
 764     void match2Dispatcher(const DevMem2D_<T>& query, const DevMem2D_<T>* trains, int n, const Mask& mask, \r
 765                           const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, \r
 766                           int cc, cudaStream_t stream)\r
 767     {\r
 768         if (query.cols <= 64)\r
 769         {\r
 770             matchUnrolledCached<16, 64, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 771         }\r
 772         else if (query.cols <= 128)\r
 773         {\r
 774             matchUnrolledCached<16, 128, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 775         }\r
 776         /*else if (query.cols <= 256)\r
 777         {\r
 778             matchUnrolled<16, 256, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 779         }\r
 780         else if (query.cols <= 512)\r
 781         {            \r
 782             matchUnrolled<16, 512, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 783         }\r
 784         else if (query.cols <= 1024)\r
 785         {            \r
 786             matchUnrolled<16, 1024, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 787         }*/\r
 788         else\r
 789         {\r
 790             match<16, Dist>(query, trains, n, mask, static_cast< DevMem2D_<int2> >(trainIdx), static_cast< DevMem2D_<int2> >(imgIdx), static_cast< DevMem2D_<float2> > (distance), stream);\r
 791         }\r
 792     }\r
 793 \r
 794     ///////////////////////////////////////////////////////////////////////////////\r
 795     // Calc distance kernel\r
 796 \r
 797     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask>\r
 798     __global__ void calcDistanceUnrolled(const DevMem2D_<T> query, const DevMem2D_<T> train, const Mask mask, PtrStepf allDist)\r
 799     {\r
 800         extern __shared__ int smem[];\r
 801 \r
 802         const int queryIdx = blockIdx.y * BLOCK_SIZE + threadIdx.y;\r
 803         const int trainIdx = blockIdx.x * BLOCK_SIZE + threadIdx.x;\r
 804 \r
 805         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 806         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 807 \r
 808         Dist dist;\r
 809 \r
 810         #pragma unroll\r
 811         for (int i = 0; i < MAX_DESC_LEN / BLOCK_SIZE; ++i)\r
 812         {\r
 813             const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 814 \r
 815             if (loadX < query.cols)\r
 816             {\r
 817                 s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = query.ptr(min(queryIdx, query.rows - 1))[loadX];\r
 818                 s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = train.ptr(min(blockIdx.x * BLOCK_SIZE + threadIdx.y, train.rows - 1))[loadX];\r
 819             }\r
 820             else\r
 821             {                \r
 822                 s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = 0;\r
 823                 s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = 0;\r
 824             }\r
 825 \r
 826             __syncthreads();\r
 827 \r
 828             #pragma unroll\r
 829             for (int j = 0; j < BLOCK_SIZE; ++j)\r
 830                 dist.reduceIter(s_query[threadIdx.y * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + threadIdx.x]);\r
 831 \r
 832             __syncthreads();\r
 833         }\r
 834 \r
 835         if (queryIdx < query.rows && trainIdx < train.rows)\r
 836         {\r
 837             float distVal = numeric_limits<float>::max();\r
 838 \r
 839             if (mask(queryIdx, trainIdx))\r
 840                 distVal = (typename Dist::result_type)dist;\r
 841 \r
 842             allDist.ptr(queryIdx)[trainIdx] = distVal;\r
 843         }\r
 844     }\r
 845 \r
 846     template <int BLOCK_SIZE, int MAX_DESC_LEN, typename Dist, typename T, typename Mask> \r
 847     void calcDistanceUnrolled(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, const DevMem2Df& allDist, cudaStream_t stream)\r
 848     {\r
 849         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 850         const dim3 grid(divUp(train.rows, BLOCK_SIZE), divUp(query.rows, BLOCK_SIZE));\r
 851 \r
 852         const size_t smemSize = (2 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 853 \r
 854         calcDistanceUnrolled<BLOCK_SIZE, MAX_DESC_LEN, Dist><<<grid, block, smemSize, stream>>>(query, train, mask, allDist);\r
 855         cudaSafeCall( cudaGetLastError() );\r
 856 \r
 857         if (stream == 0)\r
 858             cudaSafeCall( cudaDeviceSynchronize() );\r
 859     }\r
 860 \r
 861     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask>\r
 862     __global__ void calcDistance(const DevMem2D_<T> query, const DevMem2D_<T> train, const Mask mask, PtrStepf allDist)\r
 863     {\r
 864         extern __shared__ int smem[];\r
 865 \r
 866         const int queryIdx = blockIdx.y * BLOCK_SIZE + threadIdx.y;\r
 867         const int trainIdx = blockIdx.x * BLOCK_SIZE + threadIdx.x;\r
 868 \r
 869         typename Dist::value_type* s_query = (typename Dist::value_type*)(smem);\r
 870         typename Dist::value_type* s_train = (typename Dist::value_type*)(smem + BLOCK_SIZE * BLOCK_SIZE);\r
 871 \r
 872         Dist dist;\r
 873 \r
 874         for (int i = 0, endi = (query.cols + BLOCK_SIZE - 1) / BLOCK_SIZE; i < endi; ++i)\r
 875         {\r
 876             const int loadX = threadIdx.x + i * BLOCK_SIZE;\r
 877 \r
 878             if (loadX < query.cols)\r
 879             {\r
 880                 s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = query.ptr(min(queryIdx, query.rows - 1))[loadX];\r
 881                 s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = train.ptr(min(blockIdx.x * BLOCK_SIZE + threadIdx.y, train.rows - 1))[loadX];\r
 882             }\r
 883             else\r
 884             {                \r
 885                 s_query[threadIdx.y * BLOCK_SIZE + threadIdx.x] = 0;\r
 886                 s_train[threadIdx.x * BLOCK_SIZE + threadIdx.y] = 0;\r
 887             }\r
 888 \r
 889             __syncthreads();\r
 890 \r
 891             #pragma unroll\r
 892             for (int j = 0; j < BLOCK_SIZE; ++j)\r
 893                 dist.reduceIter(s_query[threadIdx.y * BLOCK_SIZE + j], s_train[j * BLOCK_SIZE + threadIdx.x]);\r
 894 \r
 895             __syncthreads();\r
 896         }\r
 897 \r
 898         if (queryIdx < query.rows && trainIdx < train.rows)\r
 899         {\r
 900             float distVal = numeric_limits<float>::max();\r
 901 \r
 902             if (mask(queryIdx, trainIdx))\r
 903                 distVal = (typename Dist::result_type)dist;\r
 904 \r
 905             allDist.ptr(queryIdx)[trainIdx] = distVal;\r
 906         }\r
 907     }\r
 908 \r
 909     template <int BLOCK_SIZE, typename Dist, typename T, typename Mask> \r
 910     void calcDistance(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, const DevMem2Df& allDist, cudaStream_t stream)\r
 911     {\r
 912         const dim3 block(BLOCK_SIZE, BLOCK_SIZE);\r
 913         const dim3 grid(divUp(train.rows, BLOCK_SIZE), divUp(query.rows, BLOCK_SIZE));\r
 914 \r
 915         const size_t smemSize = (2 * BLOCK_SIZE * BLOCK_SIZE) * sizeof(int);\r
 916 \r
 917         calcDistance<BLOCK_SIZE, Dist><<<grid, block, smemSize, stream>>>(query, train, mask, allDist);\r
 918         cudaSafeCall( cudaGetLastError() );\r
 919 \r
 920         if (stream == 0)\r
 921             cudaSafeCall( cudaDeviceSynchronize() );\r
 922     }\r
 923 \r
 924     ///////////////////////////////////////////////////////////////////////////////\r
 925     // Calc Distance dispatcher\r
 926 \r
 927     template <typename Dist, typename T, typename Mask> \r
 928     void calcDistanceDispatcher(const DevMem2D_<T>& query, const DevMem2D_<T>& train, const Mask& mask, \r
 929                                 const DevMem2Df& allDist, \r
 930                                 int cc, cudaStream_t stream)\r
 931     {\r
 932         if (query.cols <= 64)\r
 933         {\r
 934             calcDistanceUnrolled<16, 64, Dist>(query, train, mask, allDist, stream);\r
 935         }\r
 936         else if (query.cols <= 128)\r
 937         {\r
 938             calcDistanceUnrolled<16, 128, Dist>(query, train, mask, allDist, stream);\r
 939         }\r
 940         /*else if (query.cols <= 256)\r
 941         {\r
 942             calcDistanceUnrolled<16, 256, Dist>(query, train, mask, allDist, stream);\r
 943         }\r
 944         else if (query.cols <= 512)\r
 945         {            \r
 946             calcDistanceUnrolled<16, 512, Dist>(query, train, mask, allDist, stream);\r
 947         }\r
 948         else if (query.cols <= 1024)\r
 949         {            \r
 950             calcDistanceUnrolled<16, 1024, Dist>(query, train, mask, allDist, stream);\r
 951         }*/\r
 952         else\r
 953         {\r
 954             calcDistance<16, Dist>(query, train, mask, allDist, stream);\r
 955         }\r
 956     }\r
 957 \r
 958     ///////////////////////////////////////////////////////////////////////////////\r
 959     // find knn match kernel\r
 960 \r
 961     template <int BLOCK_SIZE> \r
 962     __global__ void findBestMatch(DevMem2Df allDist, int i, PtrStepi trainIdx, PtrStepf distance)\r
 963     {\r
 964         const int SMEM_SIZE = BLOCK_SIZE > 64 ? BLOCK_SIZE : 64;\r
 965         __shared__ float s_dist[SMEM_SIZE];\r
 966         __shared__ int s_trainIdx[SMEM_SIZE];\r
 967 \r
 968         const int queryIdx = blockIdx.x;\r
 969 \r
 970         float* allDistRow = allDist.ptr(queryIdx);\r
 971 \r
 972         float dist = numeric_limits<float>::max();\r
 973         int bestIdx = -1;\r
 974         \r
 975         for (int i = threadIdx.x; i < allDist.cols; i += BLOCK_SIZE)\r
 976         {\r
 977             float reg = allDistRow[i];\r
 978             if (reg < dist)\r
 979             {\r
 980                 dist = reg;\r
 981                 bestIdx = i;\r
 982             }\r
 983         }\r
 984 \r
 985         s_dist[threadIdx.x] = dist;\r
 986         s_trainIdx[threadIdx.x] = bestIdx;\r
 987         __syncthreads();\r
 988 \r
 989         reducePredVal<BLOCK_SIZE>(s_dist, dist, s_trainIdx, bestIdx, threadIdx.x, less<volatile float>());\r
 990 \r
 991         if (threadIdx.x == 0)\r
 992         {\r
 993             if (dist < numeric_limits<float>::max())\r
 994             {\r
 995                 allDistRow[bestIdx] = numeric_limits<float>::max();\r
 996                 trainIdx.ptr(queryIdx)[i] = bestIdx;\r
 997                 distance.ptr(queryIdx)[i] = dist;\r
 998             }\r
 999         }\r
1000     }\r
1001 \r
1002     template <int BLOCK_SIZE> \r
1003     void findKnnMatch(int k, const DevMem2Di& trainIdx, const DevMem2Df& distance, const DevMem2Df& allDist, cudaStream_t stream)\r
1004     {\r
1005         const dim3 block(BLOCK_SIZE, 1, 1);\r
1006         const dim3 grid(trainIdx.rows, 1, 1);\r
1007 \r
1008         for (int i = 0; i < k; ++i)\r
1009         {\r
1010             findBestMatch<BLOCK_SIZE><<<grid, block, 0, stream>>>(allDist, i, trainIdx, distance);\r
1011             cudaSafeCall( cudaGetLastError() );\r
1012         }\r
1013 \r
1014         if (stream == 0)\r
1015             cudaSafeCall( cudaDeviceSynchronize() );\r
1016     }\r
1017 \r
1018     void findKnnMatchDispatcher(int k, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream)\r
1019     {\r
1020         findKnnMatch<256>(k, static_cast<DevMem2Di>(trainIdx), static_cast<DevMem2Df>(distance), allDist, stream);\r
1021     }\r
1022 \r
1023     ///////////////////////////////////////////////////////////////////////////////\r
1024     // knn match Dispatcher\r
1025 \r
1026     template <typename Dist, typename T, typename Mask>\r
1027     void matchDispatcher(const DevMem2D_<T>& query, const DevMem2D_<T>& train, int k, const Mask& mask, \r
1028         const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, \r
1029         int cc, cudaStream_t stream)\r
1030     {\r
1031         if (k == 2)\r
1032         {\r
1033             match2Dispatcher<Dist>(query, train, mask, trainIdx, distance, cc, stream);\r
1034         }\r
1035         else\r
1036         {\r
1037             calcDistanceDispatcher<Dist>(query, train, mask, allDist, cc, stream);\r
1038             findKnnMatchDispatcher(k, trainIdx, distance, allDist, cc, stream);\r
1039         }\r
1040     }     \r
1041     \r
1042     ///////////////////////////////////////////////////////////////////////////////\r
1043     // knn match caller\r
1044 \r
1045     template <typename T> void matchL1_gpu(const DevMem2Db& query, const DevMem2Db& train, int k, const DevMem2Db& mask, \r
1046         const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, \r
1047         int cc, cudaStream_t stream)\r
1048     {\r
1049         if (mask.data)\r
1050             matchDispatcher< L1Dist<T> >(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, SingleMask(mask), trainIdx, distance, allDist, cc, stream);\r
1051         else\r
1052             matchDispatcher< L1Dist<T> >(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, WithOutMask(), trainIdx, distance, allDist, cc, stream);\r
1053     }\r
1054 \r
1055     template void matchL1_gpu<uchar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1056     //template void matchL1_gpu<schar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1057     template void matchL1_gpu<ushort>(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1058     template void matchL1_gpu<short >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1059     template void matchL1_gpu<int   >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1060     template void matchL1_gpu<float >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1061 \r
1062     template <typename T> void matchL2_gpu(const DevMem2Db& query, const DevMem2Db& train, int k, const DevMem2Db& mask, \r
1063         const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist,\r
1064         int cc, cudaStream_t stream)\r
1065     {\r
1066         if (mask.data)\r
1067             matchDispatcher<L2Dist>(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, SingleMask(mask), trainIdx, distance, allDist, cc, stream);\r
1068         else\r
1069             matchDispatcher<L2Dist>(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, WithOutMask(), trainIdx, distance, allDist, cc, stream);\r
1070     }\r
1071 \r
1072     //template void matchL2_gpu<uchar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1073     //template void matchL2_gpu<schar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1074     //template void matchL2_gpu<ushort>(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1075     //template void matchL2_gpu<short >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1076     //template void matchL2_gpu<int   >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1077     template void matchL2_gpu<float >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1078 \r
1079     template <typename T> void matchHamming_gpu(const DevMem2Db& query, const DevMem2Db& train, int k, const DevMem2Db& mask,\r
1080         const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, \r
1081         int cc, cudaStream_t stream)\r
1082     {\r
1083         if (mask.data)\r
1084             matchDispatcher<HammingDist>(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, SingleMask(mask), trainIdx, distance, allDist, cc, stream);\r
1085         else\r
1086             matchDispatcher<HammingDist>(static_cast< DevMem2D_<T> >(query), static_cast< DevMem2D_<T> >(train), k, WithOutMask(), trainIdx, distance, allDist, cc, stream);\r
1087     }\r
1088 \r
1089     template void matchHamming_gpu<uchar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1090     //template void matchHamming_gpu<schar >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1091     template void matchHamming_gpu<ushort>(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1092     //template void matchHamming_gpu<short >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1093     template void matchHamming_gpu<int   >(const DevMem2Db& queryDescs, const DevMem2Db& trainDescs, int k, const DevMem2Db& mask, const DevMem2Db& trainIdx, const DevMem2Db& distance, const DevMem2Df& allDist, int cc, cudaStream_t stream);\r
1094 \r
1095     template <typename T> void match2L1_gpu(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, \r
1096         const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, \r
1097         int cc, cudaStream_t stream)\r
1098     {\r
1099         if (masks.data)\r
1100             match2Dispatcher< L1Dist<T> >(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, MaskCollection(masks.data), trainIdx, imgIdx, distance, cc, stream);\r
1101         else\r
1102             match2Dispatcher< L1Dist<T> >(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, WithOutMask(), trainIdx, imgIdx, distance, cc, stream);\r
1103     }\r
1104     \r
1105     template void match2L1_gpu<uchar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1106     //template void match2L1_gpu<schar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1107     template void match2L1_gpu<ushort>(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1108     template void match2L1_gpu<short >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1109     template void match2L1_gpu<int   >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1110     template void match2L1_gpu<float >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1111 \r
1112     template <typename T> void match2L2_gpu(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, \r
1113         const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, \r
1114         int cc, cudaStream_t stream)\r
1115     {\r
1116         if (masks.data)\r
1117             match2Dispatcher<L2Dist>(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, MaskCollection(masks.data), trainIdx, imgIdx, distance, cc, stream);\r
1118         else\r
1119             match2Dispatcher<L2Dist>(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, WithOutMask(), trainIdx, imgIdx, distance, cc, stream);\r
1120     }\r
1121     \r
1122     //template void match2L2_gpu<uchar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1123     //template void match2L2_gpu<schar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1124     //template void match2L2_gpu<ushort>(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1125     //template void match2L2_gpu<short >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1126     //template void match2L2_gpu<int   >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Di& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1127     template void match2L2_gpu<float >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1128     \r
1129     template <typename T> void match2Hamming_gpu(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, \r
1130         const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, \r
1131         int cc, cudaStream_t stream)\r
1132     {\r
1133         if (masks.data)\r
1134             match2Dispatcher<HammingDist>(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, MaskCollection(masks.data), trainIdx, imgIdx, distance, cc, stream);\r
1135         else\r
1136             match2Dispatcher<HammingDist>(static_cast< DevMem2D_<T> >(query), (const DevMem2D_<T>*)trains.ptr(), trains.cols, WithOutMask(), trainIdx, imgIdx, distance, cc, stream);\r
1137     }\r
1138     \r
1139     template void match2Hamming_gpu<uchar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1140     //template void match2Hamming_gpu<schar >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1141     template void match2Hamming_gpu<ushort>(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1142     //template void match2Hamming_gpu<short >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1143     template void match2Hamming_gpu<int   >(const DevMem2Db& query, const DevMem2Db& trains, const DevMem2D_<PtrStepb>& masks, const DevMem2Db& trainIdx, const DevMem2Db& imgIdx, const DevMem2Db& distance, int cc, cudaStream_t stream);\r
1144 }}}\r