integral.hpp 20 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629
  1. /*M///////////////////////////////////////////////////////////////////////////////////////
  2. //
  3. // IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
  4. //
  5. // By downloading, copying, installing or using the software you agree to this license.
  6. // If you do not agree to this license, do not download, install,
  7. // copy or use the software.
  8. //
  9. //
  10. // License Agreement
  11. // For Open Source Computer Vision Library
  12. //
  13. // Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
  14. // Copyright (C) 2009, Willow Garage Inc., all rights reserved.
  15. // Copyright (C) 2013, OpenCV Foundation, all rights reserved.
  16. // Third party copyrights are property of their respective owners.
  17. //
  18. // Redistribution and use in source and binary forms, with or without modification,
  19. // are permitted provided that the following conditions are met:
  20. //
  21. // * Redistribution's of source code must retain the above copyright notice,
  22. // this list of conditions and the following disclaimer.
  23. //
  24. // * Redistribution's in binary form must reproduce the above copyright notice,
  25. // this list of conditions and the following disclaimer in the documentation
  26. // and/or other materials provided with the distribution.
  27. //
  28. // * The name of the copyright holders may not be used to endorse or promote products
  29. // derived from this software without specific prior written permission.
  30. //
  31. // This software is provided by the copyright holders and contributors "as is" and
  32. // any express or implied warranties, including, but not limited to, the implied
  33. // warranties of merchantability and fitness for a particular purpose are disclaimed.
  34. // In no event shall the Intel Corporation or contributors be liable for any direct,
  35. // indirect, incidental, special, exemplary, or consequential damages
  36. // (including, but not limited to, procurement of substitute goods or services;
  37. // loss of use, data, or profits; or business interruption) however caused
  38. // and on any theory of liability, whether in contract, strict liability,
  39. // or tort (including negligence or otherwise) arising in any way out of
  40. // the use of this software, even if advised of the possibility of such damage.
  41. //
  42. //M*/
  43. #pragma once
  44. #ifndef OPENCV_CUDEV_GRID_INTEGRAL_DETAIL_HPP
  45. #define OPENCV_CUDEV_GRID_INTEGRAL_DETAIL_HPP
  46. #include "../../common.hpp"
  47. #include "../../warp/shuffle.hpp"
  48. #include "../../block/scan.hpp"
  49. #include "../../ptr2d/glob.hpp"
  50. namespace cv { namespace cudev {
  51. namespace integral_detail
  52. {
  53. // horizontal_pass
  54. template <int NUM_SCAN_THREADS, class SrcPtr, typename D>
  55. __global__ void horizontal_pass(const SrcPtr src, GlobPtr<D> dst, const int cols)
  56. {
  57. __shared__ D smem[NUM_SCAN_THREADS * 2];
  58. __shared__ D carryElem;
  59. if (threadIdx.x == 0)
  60. carryElem = 0;
  61. __syncthreads();
  62. D* dst_row = dst.row(blockIdx.x);
  63. int numBuckets = divUp(cols, NUM_SCAN_THREADS);
  64. int offsetX = 0;
  65. while (numBuckets--)
  66. {
  67. const int curElemOffs = offsetX + threadIdx.x;
  68. D curElem = 0.0f;
  69. if (curElemOffs < cols)
  70. curElem = src(blockIdx.x, curElemOffs);
  71. const D curScanElem = blockScanInclusive<NUM_SCAN_THREADS>(curElem, smem, threadIdx.x);
  72. if (curElemOffs < cols)
  73. dst_row[curElemOffs] = carryElem + curScanElem;
  74. offsetX += NUM_SCAN_THREADS;
  75. __syncthreads();
  76. if (threadIdx.x == NUM_SCAN_THREADS - 1)
  77. {
  78. carryElem += curScanElem;
  79. }
  80. __syncthreads();
  81. }
  82. }
  83. template <int NUM_SCAN_THREADS, typename T, typename D>
  84. __global__ void horizontal_pass(const GlobPtr<T> src, GlobPtr<D> dst, const int cols)
  85. {
  86. __shared__ D smem[NUM_SCAN_THREADS * 2];
  87. __shared__ D carryElem;
  88. if (threadIdx.x == 0)
  89. carryElem = 0;
  90. __syncthreads();
  91. const T* src_row = src.row(blockIdx.x);
  92. D* dst_row = dst.row(blockIdx.x);
  93. int numBuckets = divUp(cols, NUM_SCAN_THREADS);
  94. int offsetX = 0;
  95. while (numBuckets--)
  96. {
  97. const int curElemOffs = offsetX + threadIdx.x;
  98. D curElem = 0.0f;
  99. if (curElemOffs < cols)
  100. curElem = src_row[curElemOffs];
  101. const D curScanElem = blockScanInclusive<NUM_SCAN_THREADS>(curElem, smem, threadIdx.x);
  102. if (curElemOffs < cols)
  103. dst_row[curElemOffs] = carryElem + curScanElem;
  104. offsetX += NUM_SCAN_THREADS;
  105. __syncthreads();
  106. if (threadIdx.x == NUM_SCAN_THREADS - 1)
  107. {
  108. carryElem += curScanElem;
  109. }
  110. __syncthreads();
  111. }
  112. }
  113. template <class SrcPtr, typename D>
  114. __host__ void horizontal_pass(const SrcPtr& src, const GlobPtr<D>& dst, int rows, int cols, cudaStream_t stream)
  115. {
  116. const int NUM_SCAN_THREADS = 256;
  117. const dim3 block(NUM_SCAN_THREADS);
  118. const dim3 grid(rows);
  119. horizontal_pass<NUM_SCAN_THREADS><<<grid, block, 0, stream>>>(src, dst, cols);
  120. CV_CUDEV_SAFE_CALL( cudaGetLastError() );
  121. }
  122. // horisontal_pass_8u_shfl
  123. __device__ static uchar4 int_to_uchar4(unsigned int in)
  124. {
  125. uchar4 bytes;
  126. bytes.x = (in & 0x000000ff) >> 0;
  127. bytes.y = (in & 0x0000ff00) >> 8;
  128. bytes.z = (in & 0x00ff0000) >> 16;
  129. bytes.w = (in & 0xff000000) >> 24;
  130. return bytes;
  131. }
  132. __global__ static void horisontal_pass_8u_shfl_kernel(const GlobPtr<uint4> img, GlobPtr<uint4> integral)
  133. {
  134. #if CV_CUDEV_ARCH >= 300
  135. __shared__ int sums[128];
  136. const int id = threadIdx.x;
  137. const int lane_id = id % warpSize;
  138. const int warp_id = id / warpSize;
  139. const uint4 data = img(blockIdx.x, id);
  140. const uchar4 a = int_to_uchar4(data.x);
  141. const uchar4 b = int_to_uchar4(data.y);
  142. const uchar4 c = int_to_uchar4(data.z);
  143. const uchar4 d = int_to_uchar4(data.w);
  144. int result[16];
  145. result[0] = a.x;
  146. result[1] = result[0] + a.y;
  147. result[2] = result[1] + a.z;
  148. result[3] = result[2] + a.w;
  149. result[4] = result[3] + b.x;
  150. result[5] = result[4] + b.y;
  151. result[6] = result[5] + b.z;
  152. result[7] = result[6] + b.w;
  153. result[8] = result[7] + c.x;
  154. result[9] = result[8] + c.y;
  155. result[10] = result[9] + c.z;
  156. result[11] = result[10] + c.w;
  157. result[12] = result[11] + d.x;
  158. result[13] = result[12] + d.y;
  159. result[14] = result[13] + d.z;
  160. result[15] = result[14] + d.w;
  161. int sum = result[15];
  162. // the prefix sum for each thread's 16 value is computed,
  163. // now the final sums (result[15]) need to be shared
  164. // with the other threads and add. To do this,
  165. // the shfl_up() instruction is used and a shuffle scan
  166. // operation is performed to distribute the sums to the correct
  167. // threads
  168. #pragma unroll
  169. for (int i = 1; i < 32; i *= 2)
  170. {
  171. const int n = compatible_shfl_up(sum, i, 32);
  172. if (lane_id >= i)
  173. {
  174. #pragma unroll
  175. for (int k = 0; k < 16; ++k)
  176. result[k] += n;
  177. sum += n;
  178. }
  179. }
  180. // Now the final sum for the warp must be shared
  181. // between warps. This is done by each warp
  182. // having a thread store to shared memory, then
  183. // having some other warp load the values and
  184. // compute a prefix sum, again by using shfl_up.
  185. // The results are uniformly added back to the warps.
  186. // last thread in the warp holding sum of the warp
  187. // places that in shared
  188. if (threadIdx.x % warpSize == warpSize - 1)
  189. sums[warp_id] = result[15];
  190. __syncthreads();
  191. if (warp_id == 0)
  192. {
  193. int warp_sum = sums[lane_id];
  194. #pragma unroll
  195. for (int i = 1; i < 32; i *= 2)
  196. {
  197. const int n = compatible_shfl_up(warp_sum, i, 32);
  198. if (lane_id >= i)
  199. warp_sum += n;
  200. }
  201. sums[lane_id] = warp_sum;
  202. }
  203. __syncthreads();
  204. int blockSum = 0;
  205. // fold in unused warp
  206. if (warp_id > 0)
  207. {
  208. blockSum = sums[warp_id - 1];
  209. #pragma unroll
  210. for (int k = 0; k < 16; ++k)
  211. result[k] += blockSum;
  212. }
  213. // assemble result
  214. // Each thread has 16 values to write, which are
  215. // now integer data (to avoid overflow). Instead of
  216. // each thread writing consecutive uint4s, the
  217. // approach shown here experiments using
  218. // the shuffle command to reformat the data
  219. // inside the registers so that each thread holds
  220. // consecutive data to be written so larger contiguous
  221. // segments can be assembled for writing.
  222. /*
  223. For example data that needs to be written as
  224. GMEM[16] <- x0 x1 x2 x3 y0 y1 y2 y3 z0 z1 z2 z3 w0 w1 w2 w3
  225. but is stored in registers (r0..r3), in four threads (0..3) as:
  226. threadId 0 1 2 3
  227. r0 x0 y0 z0 w0
  228. r1 x1 y1 z1 w1
  229. r2 x2 y2 z2 w2
  230. r3 x3 y3 z3 w3
  231. after apply shfl_xor operations to move data between registers r1..r3:
  232. threadId 00 01 10 11
  233. x0 y0 z0 w0
  234. xor(01)->y1 x1 w1 z1
  235. xor(10)->z2 w2 x2 y2
  236. xor(11)->w3 z3 y3 x3
  237. and now x0..x3, and z0..z3 can be written out in order by all threads.
  238. In the current code, each register above is actually representing
  239. four integers to be written as uint4's to GMEM.
  240. */
  241. result[4] = shfl_xor(result[4] , 1, 32);
  242. result[5] = shfl_xor(result[5] , 1, 32);
  243. result[6] = shfl_xor(result[6] , 1, 32);
  244. result[7] = shfl_xor(result[7] , 1, 32);
  245. result[8] = shfl_xor(result[8] , 2, 32);
  246. result[9] = shfl_xor(result[9] , 2, 32);
  247. result[10] = shfl_xor(result[10], 2, 32);
  248. result[11] = shfl_xor(result[11], 2, 32);
  249. result[12] = shfl_xor(result[12], 3, 32);
  250. result[13] = shfl_xor(result[13], 3, 32);
  251. result[14] = shfl_xor(result[14], 3, 32);
  252. result[15] = shfl_xor(result[15], 3, 32);
  253. uint4* integral_row = integral.row(blockIdx.x);
  254. uint4 output;
  255. ///////
  256. if (threadIdx.x % 4 == 0)
  257. output = make_uint4(result[0], result[1], result[2], result[3]);
  258. if (threadIdx.x % 4 == 1)
  259. output = make_uint4(result[4], result[5], result[6], result[7]);
  260. if (threadIdx.x % 4 == 2)
  261. output = make_uint4(result[8], result[9], result[10], result[11]);
  262. if (threadIdx.x % 4 == 3)
  263. output = make_uint4(result[12], result[13], result[14], result[15]);
  264. integral_row[threadIdx.x % 4 + (threadIdx.x / 4) * 16] = output;
  265. ///////
  266. if (threadIdx.x % 4 == 2)
  267. output = make_uint4(result[0], result[1], result[2], result[3]);
  268. if (threadIdx.x % 4 == 3)
  269. output = make_uint4(result[4], result[5], result[6], result[7]);
  270. if (threadIdx.x % 4 == 0)
  271. output = make_uint4(result[8], result[9], result[10], result[11]);
  272. if (threadIdx.x % 4 == 1)
  273. output = make_uint4(result[12], result[13], result[14], result[15]);
  274. integral_row[(threadIdx.x + 2) % 4 + (threadIdx.x / 4) * 16 + 8] = output;
  275. // continuning from the above example,
  276. // this use of shfl_xor() places the y0..y3 and w0..w3 data
  277. // in order.
  278. #pragma unroll
  279. for (int i = 0; i < 16; ++i)
  280. result[i] = shfl_xor(result[i], 1, 32);
  281. if (threadIdx.x % 4 == 0)
  282. output = make_uint4(result[0], result[1], result[2], result[3]);
  283. if (threadIdx.x % 4 == 1)
  284. output = make_uint4(result[4], result[5], result[6], result[7]);
  285. if (threadIdx.x % 4 == 2)
  286. output = make_uint4(result[8], result[9], result[10], result[11]);
  287. if (threadIdx.x % 4 == 3)
  288. output = make_uint4(result[12], result[13], result[14], result[15]);
  289. integral_row[threadIdx.x % 4 + (threadIdx.x / 4) * 16 + 4] = output;
  290. ///////
  291. if (threadIdx.x % 4 == 2)
  292. output = make_uint4(result[0], result[1], result[2], result[3]);
  293. if (threadIdx.x % 4 == 3)
  294. output = make_uint4(result[4], result[5], result[6], result[7]);
  295. if (threadIdx.x % 4 == 0)
  296. output = make_uint4(result[8], result[9], result[10], result[11]);
  297. if (threadIdx.x % 4 == 1)
  298. output = make_uint4(result[12], result[13], result[14], result[15]);
  299. integral_row[(threadIdx.x + 2) % 4 + (threadIdx.x / 4) * 16 + 12] = output;
  300. #endif
  301. }
  302. __host__ static void horisontal_pass_8u_shfl(const GlobPtr<uchar> src, GlobPtr<uint> integral, int rows, int cols, cudaStream_t stream)
  303. {
  304. // each thread handles 16 values, use 1 block/row
  305. // save, because step is actually can't be less 512 bytes
  306. const int block = cols / 16;
  307. // launch 1 block / row
  308. const int grid = rows;
  309. CV_CUDEV_SAFE_CALL( cudaFuncSetCacheConfig(horisontal_pass_8u_shfl_kernel, cudaFuncCachePreferL1) );
  310. GlobPtr<uint4> src4 = globPtr((uint4*) src.data, src.step);
  311. GlobPtr<uint4> integral4 = globPtr((uint4*) integral.data, integral.step);
  312. horisontal_pass_8u_shfl_kernel<<<grid, block, 0, stream>>>(src4, integral4);
  313. CV_CUDEV_SAFE_CALL( cudaGetLastError() );
  314. }
  315. // vertical
  316. template <typename T>
  317. __global__ void vertical_pass(GlobPtr<T> integral, const int rows, const int cols)
  318. {
  319. #if CV_CUDEV_ARCH >= 300
  320. __shared__ T sums[32][9];
  321. const int tidx = blockIdx.x * blockDim.x + threadIdx.x;
  322. const int lane_id = tidx % 8;
  323. sums[threadIdx.x][threadIdx.y] = 0;
  324. __syncthreads();
  325. T stepSum = 0;
  326. int numBuckets = divUp(rows, blockDim.y);
  327. int y = threadIdx.y;
  328. while (numBuckets--)
  329. {
  330. T* p = integral.row(y) + tidx;
  331. T sum = (tidx < cols) && (y < rows) ? *p : 0;
  332. sums[threadIdx.x][threadIdx.y] = sum;
  333. __syncthreads();
  334. // place into SMEM
  335. // shfl scan reduce the SMEM, reformating so the column
  336. // sums are computed in a warp
  337. // then read out properly
  338. const int j = threadIdx.x % 8;
  339. const int k = threadIdx.x / 8 + threadIdx.y * 4;
  340. T partial_sum = sums[k][j];
  341. for (int i = 1; i <= 8; i *= 2)
  342. {
  343. T n = compatible_shfl_up(partial_sum, i, 32);
  344. if (lane_id >= i)
  345. partial_sum += n;
  346. }
  347. sums[k][j] = partial_sum;
  348. __syncthreads();
  349. if (threadIdx.y > 0)
  350. sum += sums[threadIdx.x][threadIdx.y - 1];
  351. sum += stepSum;
  352. stepSum += sums[threadIdx.x][blockDim.y - 1];
  353. __syncthreads();
  354. if ((tidx < cols) && (y < rows))
  355. {
  356. *p = sum;
  357. }
  358. y += blockDim.y;
  359. }
  360. #else
  361. __shared__ T smem[32][32];
  362. __shared__ T prevVals[32];
  363. volatile T* smem_row = &smem[0][0] + 64 * threadIdx.y;
  364. if (threadIdx.y == 0)
  365. prevVals[threadIdx.x] = 0;
  366. __syncthreads();
  367. const int x = blockIdx.x * blockDim.x + threadIdx.x;
  368. int numBuckets = divUp(rows, 8 * 4);
  369. int offsetY = 0;
  370. while (numBuckets--)
  371. {
  372. const int curRowOffs = offsetY + threadIdx.y;
  373. T curElems[4];
  374. T temp[4];
  375. // load patch
  376. smem[threadIdx.y + 0][threadIdx.x] = 0.0f;
  377. smem[threadIdx.y + 8][threadIdx.x] = 0.0f;
  378. smem[threadIdx.y + 16][threadIdx.x] = 0.0f;
  379. smem[threadIdx.y + 24][threadIdx.x] = 0.0f;
  380. if (x < cols)
  381. {
  382. for (int i = 0; i < 4; ++i)
  383. {
  384. if (curRowOffs + i * 8 < rows)
  385. smem[threadIdx.y + i * 8][threadIdx.x] = integral(curRowOffs + i * 8, x);
  386. }
  387. }
  388. __syncthreads();
  389. // reduce
  390. curElems[0] = smem[threadIdx.x][threadIdx.y ];
  391. curElems[1] = smem[threadIdx.x][threadIdx.y + 8];
  392. curElems[2] = smem[threadIdx.x][threadIdx.y + 16];
  393. curElems[3] = smem[threadIdx.x][threadIdx.y + 24];
  394. __syncthreads();
  395. temp[0] = curElems[0] = warpScanInclusive(curElems[0], smem_row, threadIdx.x);
  396. temp[1] = curElems[1] = warpScanInclusive(curElems[1], smem_row, threadIdx.x);
  397. temp[2] = curElems[2] = warpScanInclusive(curElems[2], smem_row, threadIdx.x);
  398. temp[3] = curElems[3] = warpScanInclusive(curElems[3], smem_row, threadIdx.x);
  399. curElems[0] += prevVals[threadIdx.y ];
  400. curElems[1] += prevVals[threadIdx.y + 8];
  401. curElems[2] += prevVals[threadIdx.y + 16];
  402. curElems[3] += prevVals[threadIdx.y + 24];
  403. __syncthreads();
  404. if (threadIdx.x == 31)
  405. {
  406. prevVals[threadIdx.y ] += temp[0];
  407. prevVals[threadIdx.y + 8] += temp[1];
  408. prevVals[threadIdx.y + 16] += temp[2];
  409. prevVals[threadIdx.y + 24] += temp[3];
  410. }
  411. smem[threadIdx.y ][threadIdx.x] = curElems[0];
  412. smem[threadIdx.y + 8][threadIdx.x] = curElems[1];
  413. smem[threadIdx.y + 16][threadIdx.x] = curElems[2];
  414. smem[threadIdx.y + 24][threadIdx.x] = curElems[3];
  415. __syncthreads();
  416. // store patch
  417. if (x < cols)
  418. {
  419. // read 4 value from source
  420. for (int i = 0; i < 4; ++i)
  421. {
  422. if (curRowOffs + i * 8 < rows)
  423. integral(curRowOffs + i * 8, x) = smem[threadIdx.x][threadIdx.y + i * 8];
  424. }
  425. }
  426. __syncthreads();
  427. offsetY += 8 * 4;
  428. }
  429. #endif
  430. }
  431. template <typename T>
  432. __host__ void vertical_pass(const GlobPtr<T>& integral, int rows, int cols, cudaStream_t stream)
  433. {
  434. const dim3 block(32, 8);
  435. const dim3 grid(divUp(cols, block.x));
  436. vertical_pass<<<grid, block, 0, stream>>>(integral, rows, cols);
  437. CV_CUDEV_SAFE_CALL( cudaGetLastError() );
  438. }
  439. // integral
  440. template <class SrcPtr, typename D>
  441. __host__ void integral(const SrcPtr& src, const GlobPtr<D>& dst, int rows, int cols, cudaStream_t stream)
  442. {
  443. horizontal_pass(src, dst, rows, cols, stream);
  444. vertical_pass(dst, rows, cols, stream);
  445. if (stream == 0)
  446. CV_CUDEV_SAFE_CALL( cudaDeviceSynchronize() );
  447. }
  448. __host__ static void integral(const GlobPtr<uchar>& src, const GlobPtr<uint>& dst, int rows, int cols, cudaStream_t stream)
  449. {
  450. if (deviceSupports(FEATURE_SET_COMPUTE_30)
  451. && (cols % 64 == 0)
  452. && reinterpret_cast<intptr_t>(src.data) % 32 == 0
  453. && reinterpret_cast<intptr_t>(dst.data) % 32 == 0)
  454. {
  455. horisontal_pass_8u_shfl(src, dst, rows, cols, stream);
  456. }
  457. else
  458. {
  459. horizontal_pass(src, dst, rows, cols, stream);
  460. }
  461. vertical_pass(dst, rows, cols, stream);
  462. if (stream == 0)
  463. CV_CUDEV_SAFE_CALL( cudaDeviceSynchronize() );
  464. }
  465. __host__ __forceinline__ void integral(const GlobPtr<uchar>& src, const GlobPtr<int>& dst, int rows, int cols, cudaStream_t stream)
  466. {
  467. GlobPtr<uint> dstui = globPtr((uint*) dst.data, dst.step);
  468. integral(src, dstui, rows, cols, stream);
  469. }
  470. }
  471. }}
  472. #endif