kernel.cu 27 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576577578579580581582583584585586587588589590591592593594595596597598599600601602603604605606607608609610611612613614615616617618619620621622623624625626627628629630631632633634635636637638639640641642643644645646647648649650651652653654655656657658659660661662663664665666667668669670671672673674675676677678679680681682683684685686687688689690691692693694695696697698699700701702703704705706707708709710711712713714715716717718719720721722723724725726727728729730731732733734735736737738739740741742743744745746747748749750751752753754755756757758759760761762763764765766767768769770771772773774775776777778779780781782783784785786787788789790791792793794795796797798799800801802803804805806807808809810811812813814815816817818819820821822823824
  1. /*
  2. %% Author: Rodrigo de Barros Vimieiro
  3. % Date: March, 2019
  4. % rodrigo.vimieiro@gmail.com
  5. % =========================================================================
  6. %{
  7. % -------------------------------------------------------------------------
  8. % projectionDDb_mex_CUDA(data3d,param)
  9. % -------------------------------------------------------------------------
  10. % DESCRIPTION:
  11. % This function calculates the volume projection based on the
  12. % Branchless Distance-Driven principle.
  13. % The geometry is for DBT with half cone-beam. All parameters are set
  14. % in "ParameterSettings" code.
  15. %
  16. % INPUT:
  17. %
  18. % - data3d = 3D volume for projection
  19. % - param = Parameter of all geometry
  20. % - nProj = projection number to be projected
  21. %
  22. % OUTPUT:
  23. %
  24. % - proj = projections for each angle.
  25. %
  26. % Reference:
  27. % - Branchless Distance Driven Projection and Backprojection,
  28. % Samit Basu and Bruno De Man (2006)
  29. % - GPU Acceleration of Branchless Distance Driven Projection and
  30. % Backprojection, Liu et al (2016)
  31. % - GPU-Based Branchless Distance-Driven Projection and Backprojection,
  32. % Liu et al (2017)
  33. % - A GPU Implementation of Distance-Driven Computed Tomography,
  34. % Ryan D. Wagner (2017)
  35. % ---------------------------------------------------------------------
  36. % Copyright (C) <2019> <Rodrigo de Barros Vimieiro>
  37. %
  38. % This program is free software: you can redistribute it and/or modify
  39. % it under the terms of the GNU General Public License as published by
  40. % the Free Software Foundation, either version 3 of the License, or
  41. % (at your option) any later version.
  42. %
  43. % This program is distributed in the hope that it will be useful,
  44. % but WITHOUT ANY WARRANTY; without even the implied warranty of
  45. % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  46. % GNU General Public License for more details.
  47. %
  48. % You should have received a copy of the GNU General Public License
  49. % along with this program. If not, see <http://www.gnu.org/licenses/>.
  50. %}
  51. % =========================================================================
  52. %% 3-D Projection Branchless Distance Driven Code
  53. */
  54. #include "projectionDDb_mex_cuda.h"
  55. /****************************************************************************
  56. * CUDA Kernels *
  57. ****************************************************************************/
  58. __global__ void pad_volume_kernel(double* d_img, const int nPixXMap, const int nPixYMap, const int nElem, unsigned int nz) {
  59. const int threadGId = blockIdx.x * blockDim.x + threadIdx.x;
  60. // Make sure we don't try and access memory outside
  61. // by having any threads mapped there return early
  62. if (threadGId >= nElem)
  63. return;
  64. d_img[(nz*nPixYMap *nPixXMap) + (threadGId*nPixYMap)] = 0;
  65. return;
  66. }
  67. __global__ void map_boudaries_kernel(double* d_pBound, const int nElem, const double valueLeftBound, const double sizeElem, const double offset) {
  68. const int threadGId = blockIdx.x * blockDim.x + threadIdx.x;
  69. // Make sure we don't try and access memory outside
  70. // by having any threads mapped there return early
  71. if (threadGId >= nElem)
  72. return;
  73. d_pBound[threadGId] = (threadGId - valueLeftBound) * sizeElem + offset;
  74. return;
  75. }
  76. __global__ void rot_detector_kernel(double* d_pRdetY, double* d_pRdetZ, double* d_pYcoord, double* d_pZcoord, const double yOffset, const double zOffset,
  77. const double phi, const int nElem) {
  78. const int threadGId = blockIdx.x * blockDim.x + threadIdx.x;
  79. // Make sure we don't try and access memory outside
  80. // by having any threads mapped there return early
  81. if (threadGId >= nElem)
  82. return;
  83. // cos and sin are in measured in radians.
  84. d_pRdetY[threadGId] = ((d_pYcoord[threadGId] - yOffset)* cos(phi) - (d_pZcoord[threadGId] - zOffset)* sin(phi)) + yOffset;
  85. d_pRdetZ[threadGId] = ((d_pYcoord[threadGId] - yOffset)* sin(phi) + (d_pZcoord[threadGId] - zOffset)* cos(phi)) + zOffset;
  86. return;
  87. }
  88. __global__ void mapDet2Slice_kernel(double* const pXmapp, double* const pYmapp, double tubeX, double tubeY, double tubeZ, double* const pXcoord,
  89. double* const pYcoord, double* const pZcoord, double* const pZSlicecoord, const int nDetXMap, const int nDetYMap, const unsigned int nz) {
  90. /*
  91. Note: Matlab has linear indexing as a column of elements, i.e, the elements are actually stored in memory as queued columns.
  92. So threads in X are launched in the column direction (DBT Y coord), and threads in Y are launched in the row direction (DBT X coord).
  93. */
  94. const int2 thread_2D_pos = make_int2(blockIdx.x * blockDim.x + threadIdx.x,
  95. blockIdx.y * blockDim.y + threadIdx.y);
  96. const int thread_1D_pos = thread_2D_pos.y * nDetYMap + thread_2D_pos.x;
  97. // Make sure we don't try and access memory outside the detector
  98. // by having any threads mapped there return early
  99. if (thread_2D_pos.x >= nDetYMap || thread_2D_pos.y >= nDetXMap)
  100. return;
  101. pXmapp[thread_1D_pos] = ((pXcoord[thread_2D_pos.y] - tubeX)*(pZSlicecoord[nz] - pZcoord[thread_2D_pos.x]) - (pXcoord[thread_2D_pos.y]*tubeZ) + (pXcoord[thread_2D_pos.y]*pZcoord[thread_2D_pos.x]))/(-tubeZ + pZcoord[thread_2D_pos.x]);
  102. if (thread_2D_pos.y == 0)
  103. pYmapp[thread_2D_pos.x] = ((pYcoord[thread_2D_pos.x] - tubeY)*(pZSlicecoord[nz] - pZcoord[thread_2D_pos.x]) - (pYcoord[thread_2D_pos.x]*tubeZ) + (pYcoord[thread_2D_pos.x]*pZcoord[thread_2D_pos.x]))/(-tubeZ + pZcoord[thread_2D_pos.x]);
  104. return;
  105. }
  106. __global__ void img_integration_kernel(double* d_img, const int nPixX, const int nPixY, bool direction, unsigned int offsetX, unsigned int offsetY, unsigned int nSlices) {
  107. /*
  108. Integration of 2D slices over the whole volume
  109. (S.1.Integration. - Liu et al(2017))
  110. ** Perfom an inclusive scan **
  111. */
  112. const int3 memory_2D_pos = make_int3(blockIdx.x * blockDim.x + threadIdx.x + offsetX,
  113. blockIdx.y * blockDim.y + threadIdx.y + offsetY,
  114. blockIdx.z * blockDim.z + threadIdx.z);
  115. const int2 thread_2D_pos = make_int2(blockIdx.x * blockDim.x + threadIdx.x,
  116. blockIdx.y * blockDim.y + threadIdx.y);
  117. // Make sure we don't try and access memory outside the detector
  118. // by having any threads mapped there return early
  119. if (memory_2D_pos.x >= nPixY || memory_2D_pos.y >= nPixX || memory_2D_pos.z >= nSlices)
  120. return;
  121. if (direction == integrateXcoord){
  122. for (int s = 1; s <= blockDim.y; s *= 2) {
  123. int spot = thread_2D_pos.y - s;
  124. double val = 0;
  125. if (spot >= 0) {
  126. val = d_img[(memory_2D_pos.z*nPixY*nPixX) + (offsetY + spot) * nPixY + memory_2D_pos.x];
  127. }
  128. __syncthreads();
  129. if (spot >= 0) {
  130. d_img[(memory_2D_pos.z*nPixY*nPixX) + (memory_2D_pos.y * nPixY) + memory_2D_pos.x] += val;
  131. }
  132. __syncthreads();
  133. }
  134. }
  135. else
  136. {
  137. for (int s = 1; s <= blockDim.x; s *= 2) {
  138. int spot = thread_2D_pos.x - s;
  139. double val = 0;
  140. if (spot >= 0) {
  141. val = d_img[(memory_2D_pos.z*nPixY*nPixX) + memory_2D_pos.y * nPixY + spot + offsetX];
  142. }
  143. __syncthreads();
  144. if (spot >= 0) {
  145. d_img[(memory_2D_pos.z*nPixY*nPixX) + (memory_2D_pos.y * nPixY) + memory_2D_pos.x] += val;
  146. }
  147. __syncthreads();
  148. }
  149. }
  150. return;
  151. }
  152. //__global__ void bilinear_interp_kernel(double* d_interpData, cudaTextureObject_t texObj, double* const d_pDetmX, double* const d_pDetmY, const int nDetXMap, const int nDetYMap, const double2 normFactor, const double2 offset)
  153. //{
  154. //
  155. // const int2 thread_2D_pos = make_int2(blockIdx.x * blockDim.x + threadIdx.x,
  156. // blockIdx.y * blockDim.y + threadIdx.y);
  157. //
  158. // // Make sure we don't try and access memory outside the detector
  159. // // by having any threads mapped there return early
  160. // if (thread_2D_pos.x >= nDetYMap || thread_2D_pos.y >= nDetXMap)
  161. // return;
  162. //
  163. // // Normalize coordinates from 0-1 (coords on texture memory)
  164. // double xQuery = d_pDetmX[thread_2D_pos.y * nDetYMap + thread_2D_pos.x] / normFactor.x + offset.x + (-.5f * (1 / normFactor.x));
  165. // double yQuery = d_pDetmY[thread_2D_pos.x] / normFactor.y + offset.y + (.5f * (1 / normFactor.y));
  166. //
  167. // // Read from texture and write to global memory
  168. // d_interpData[thread_2D_pos.y * nDetYMap + thread_2D_pos.x] = tex2D<double>(texObj, xQuery, yQuery);
  169. //
  170. // return;
  171. //}
  172. __global__ void bilinear_interpolation_kernel_GPU(double* d_projI, double* d_pVolume, double* d_pDetmX, double* d_pDetmY, const int nDetXMap, const int nDetYMap,
  173. const int nPixXMap, const int nPixYMap, const int nPixX, const int nPixY, const double pixelSize, const unsigned int nz) {
  174. const int2 thread_2D_pos = make_int2(blockIdx.x * blockDim.x + threadIdx.x,
  175. blockIdx.y * blockDim.y + threadIdx.y);
  176. // Make sure we don't try and access memory outside the detector
  177. // by having any threads mapped there return early
  178. if (thread_2D_pos.x >= nDetYMap || thread_2D_pos.y >= nDetXMap)
  179. return;
  180. /*
  181. S.2. Interpolation - Liu et al (2017)
  182. Reference:
  183. - https://en.wikipedia.org/wiki/Bilinear_interpolation
  184. - https://stackoverflow.com/questions/21128731/bilinear-interpolation-in-c-c-and-cuda
  185. Note: *** We are using the Unit Square equation ***
  186. alpha = X - X1
  187. 1-alpha = X2 - X
  188. beta = Y - Y1
  189. 1-beta = Y2 - Y
  190. ----> Xcoord (thread_2D_pos.y)
  191. | 0___________
  192. v |_d00_|_d10_|
  193. Ycoord (thread_2D_pos.x) |_d01_|_d11_|
  194. */
  195. // Adjust the mapped coordinates to cross the range of (0-nPixX).*dx
  196. // Divide by pixelSize to get a unitary pixel size
  197. const double xNormData = nPixX - d_pDetmX[thread_2D_pos.y * nDetYMap + thread_2D_pos.x] / pixelSize;
  198. const signed int xData = floor(xNormData);
  199. const double alpha = xNormData - xData;
  200. // Adjust the mapped coordinates to cross the range of (0-nPixY).*dy
  201. // Divide by pixelSize to get a unitary pixel size
  202. const double yNormData = (nPixY / 2.0) + (d_pDetmY[thread_2D_pos.x] / pixelSize);
  203. const signed int yData = floor(yNormData);
  204. const double beta = yNormData - yData;
  205. double d00, d01, d10, d11;
  206. if (((xNormData) >= 0) && ((xNormData) <= nPixX) && ((yNormData) >= 0) && ((yNormData) <= nPixY)) d00 = d_pVolume[(nz*nPixYMap*nPixXMap) + (xData*nPixYMap + yData)]; else d00 = 0.0;
  207. if (((xData + 1) > 0) && ((xData + 1) <= nPixX) && ((yNormData) >= 0) && ((yNormData) <= nPixY)) d10 = d_pVolume[(nz*nPixYMap*nPixXMap) + ((xData + 1)*nPixYMap + yData)]; else d10 = 0.0;
  208. if (((xNormData) >= 0) && ((xNormData) <= nPixX) && ((yData + 1) > 0) && ((yData + 1) <= nPixY)) d01 = d_pVolume[(nz*nPixYMap*nPixXMap) + (xData*nPixYMap + yData + 1)]; else d01 = 0.0;
  209. if (((xData + 1) > 0) && ((xData + 1) <= nPixX) && ((yData + 1) > 0) && ((yData + 1) <= nPixY)) d11 = d_pVolume[(nz*nPixYMap*nPixXMap) + ((xData + 1)*nPixYMap + yData + 1)]; else d11 = 0.0;
  210. double result_temp1 = alpha * d10 + (-d00 * alpha + d00);
  211. double result_temp2 = alpha * d11 + (-d01 * alpha + d01);
  212. d_projI[thread_2D_pos.y * nDetYMap + thread_2D_pos.x] = beta * result_temp2 + (-result_temp1 * beta + result_temp1);
  213. }
  214. __global__ void differentiation_kernel(double* d_pProj, double* d_projI, double* const d_pDetmX, double* const d_pDetmY, double tubeX, double rtubeY, double rtubeZ,
  215. double* const d_pDetX, double* const d_pRdetY, double* const d_pRdetZ, const int nDetX, const int nDetY, const int nDetXMap,
  216. const int nDetYMap, const double du, const double dv, const double dx, const double dy, const double dz, const unsigned int p) {
  217. const int2 thread_2D_pos = make_int2(blockIdx.x * blockDim.x + threadIdx.x,
  218. blockIdx.y * blockDim.y + threadIdx.y);
  219. const int thread_1D_pos = (nDetX*nDetY*p) + (thread_2D_pos.y * nDetY) + thread_2D_pos.x;
  220. // Make sure we don't try and access memory outside the detector
  221. // by having any threads mapped there return early
  222. if (thread_2D_pos.x >= nDetY || thread_2D_pos.y >= nDetX)
  223. return;
  224. /*
  225. S.3. Differentiation - Eq. 24 - Liu et al (2017)
  226. Detector integral projection
  227. ___________
  228. |_A_|_B_|___|
  229. |_C_|_D_|___|
  230. |___|___|___|
  231. (thread_2D_pos.x,thread_2D_pos.y)
  232. ________________
  233. |_A_|__B__|_____|
  234. |_C_|(0,0)|(0,1)|
  235. |___|(1,0)|(1,1)|
  236. Threads are lauched from D up to nDetX (thread_2D_pos.y) and nDetY (thread_2D_pos.x)
  237. i.e., they are running on the detector image. Thread (0,0) is on D.
  238. Coordinates on intergal projection:
  239. A = thread_2D_pos.y * nDetYMap + thread_2D_pos.x
  240. B = ((thread_2D_pos.y+1) * nDetYMap) + thread_2D_pos.x
  241. C = thread_2D_pos.y * nDetYMap + thread_2D_pos.x + 1
  242. D = ((thread_2D_pos.y+1) * nDetYMap) + thread_2D_pos.x + 1
  243. */
  244. unsigned int coordA = thread_2D_pos.y * nDetYMap + thread_2D_pos.x;
  245. unsigned int coordB = ((thread_2D_pos.y + 1) * nDetYMap) + thread_2D_pos.x;
  246. unsigned int coordC = coordA + 1;
  247. unsigned int coordD = coordB + 1;
  248. // Detector X coord overlap calculation
  249. double detMapSizeX = d_pDetmX[coordA] - d_pDetmX[coordB];
  250. // Detector Y coord overlap calculation
  251. double detMapSizeY = d_pDetmY[thread_2D_pos.x+1] - d_pDetmY[thread_2D_pos.x];
  252. // x - ray angle in X coord
  253. double gamma = atan((d_pDetX[thread_2D_pos.y] + (du / 2.0) - tubeX) / rtubeZ);
  254. // x - ray angle in Y coord
  255. double alpha = atan((d_pRdetY[thread_2D_pos.x] + (dv / 2.0) - rtubeY) / rtubeZ);
  256. double dA, dB, dC, dD;
  257. dA = d_projI[coordA];
  258. dB = d_projI[coordB];
  259. dC = d_projI[coordC];
  260. dD = d_projI[coordD];
  261. // Treat border of interpolated integral detector
  262. if (dC == 0 && dD == 0) {
  263. dC = dA;
  264. dD = dB;
  265. }
  266. // S.3.Differentiation - Eq. 24 - Liu et al(2017)
  267. d_pProj[thread_1D_pos] += ((dD - dC - dB + dA)*(dx*dy*dz / (cos(alpha)*cos(gamma)*detMapSizeX*detMapSizeY)));
  268. return;
  269. }
  270. /****************************************************************************
  271. * function: projectionDDb() - CUDA projection Branchless Distance Driven. *
  272. ****************************************************************************/
  273. void projectionDDb(double* const h_pProj,
  274. double* const h_pVolume,
  275. double* const h_pTubeAngle,
  276. double* const h_pDetAngle,
  277. const double idXProj,
  278. const unsigned int nProj,
  279. const unsigned int nPixX,
  280. const unsigned int nPixY,
  281. const unsigned int nSlices,
  282. const unsigned int nDetX,
  283. const unsigned int nDetY,
  284. const double dx,
  285. const double dy,
  286. const double dz,
  287. const double du,
  288. const double dv,
  289. const double DSD,
  290. const double DDR,
  291. const double DAG)
  292. {
  293. // Unique GPU
  294. int devID = 0;
  295. cudaDeviceProp deviceProp;
  296. cudaGetDeviceProperties(&deviceProp, devID);
  297. cudaCheckErrors("cudaGetDeviceProperties");
  298. const unsigned int maxThreadsPerBlock = deviceProp.maxThreadsPerBlock;
  299. mexPrintf("GPU Device %d: \"%s\" with compute capability %d.%d has %d Multi-Processors and %zu bytes of global memory\n\n", devID,
  300. deviceProp.name, deviceProp.major, deviceProp.minor, deviceProp.multiProcessorCount, deviceProp.totalGlobalMem);
  301. // Create timmer variables and start it
  302. StopWatchInterface *timer = NULL;
  303. sdkCreateTimer(&timer);
  304. sdkStartTimer(&timer);
  305. //cudaStream_t stream1;
  306. //cudaStreamCreate(&stream1);
  307. dim3 threadsPerBlock(1, 1, 1);
  308. dim3 blockSize(1, 1, 1);
  309. // Number of mapped detector and pixels
  310. const int nDetXMap = nDetX + 1;
  311. const int nDetYMap = nDetY + 1;
  312. const int nPixXMap = nPixX + 1;
  313. const int nPixYMap = nPixY + 1;
  314. // Convention: h_ variables live on host
  315. // Convention: d_ variables live on device (GPU global mem)
  316. double* d_pProj;
  317. double* d_projI;
  318. double* d_pVolume;
  319. double* d_pTubeAngle;
  320. double* d_pDetAngle;
  321. // Allocate global memory on the device, place result in "d_----"
  322. cudaMalloc((void **)&d_pProj, nDetX*nDetY*nProj * sizeof(double));
  323. cudaMalloc((void **)&d_projI, nDetXMap*nDetYMap * sizeof(double));
  324. cudaMalloc((void **)&d_pVolume, nPixXMap*nPixYMap*nSlices * sizeof(double));
  325. cudaMalloc((void **)&d_pTubeAngle, nProj * sizeof(double));
  326. cudaMalloc((void **)&d_pDetAngle, nProj * sizeof(double));
  327. cudaCheckErrors("cudaMalloc Initial");
  328. // Copy data from host memory "h_----" to device memory "d_----"
  329. cudaMemcpy((void *)d_pTubeAngle, (void *)h_pTubeAngle, nProj * sizeof(double), cudaMemcpyHostToDevice);
  330. cudaMemcpy((void *)d_pDetAngle, (void *)h_pDetAngle, nProj * sizeof(double), cudaMemcpyHostToDevice);
  331. cudaCheckErrors("cudaMemcpy Angles");
  332. /*
  333. Copy volume data from host memory "h_----" to device memory "d_----" padding with zeros for image integation
  334. */
  335. // Initialize first column and row with zeros
  336. double* h_pVolume_tmp;
  337. double* d_pVolume_tmp;
  338. threadsPerBlock.x = maxThreadsPerBlock;
  339. blockSize.x = (nPixXMap / maxThreadsPerBlock) + 1;
  340. for (int nz = 0; nz < nSlices; nz++) {
  341. // Pad on X coord direction
  342. pad_volume_kernel << <blockSize, threadsPerBlock >> > (d_pVolume, nPixXMap, nPixYMap, nPixXMap, nz);
  343. // Pad on Y coord direction
  344. d_pVolume_tmp = d_pVolume + (nPixXMap*nPixYMap*nz) + 1;
  345. cudaMemset(d_pVolume_tmp, 0, nPixY * sizeof(double));
  346. cudaCheckErrors("cudaMemset Padding Volume");
  347. }
  348. // Copy volume data from host memory
  349. for (unsigned int nz = 0; nz < nSlices; nz++)
  350. for (unsigned int c = 0; c < nPixX; c++) {
  351. h_pVolume_tmp = h_pVolume + (c *nPixY) + (nPixX*nPixY*nz);
  352. d_pVolume_tmp = d_pVolume + (((c+1) *nPixYMap) + 1) + (nPixXMap*nPixYMap*nz);
  353. cudaMemcpy((void *)d_pVolume_tmp, (void *)h_pVolume_tmp, nPixY * sizeof(double), cudaMemcpyHostToDevice);
  354. cudaCheckErrors("cudaMemset Volume");
  355. }
  356. // Pointer for projections coordinates
  357. double* d_pDetX;
  358. double* d_pDetY;
  359. double* d_pDetZ;
  360. double* d_pObjX;
  361. double* d_pObjY;
  362. double* d_pObjZ;
  363. // Allocate global memory on the device for projections coordinates
  364. cudaMalloc((void **)&d_pDetX, nDetXMap * sizeof(double));
  365. cudaMalloc((void **)&d_pDetY, nDetYMap * sizeof(double));
  366. cudaMalloc((void **)&d_pDetZ, nDetYMap * sizeof(double));
  367. cudaMalloc((void **)&d_pObjX, nPixXMap * sizeof(double));
  368. cudaMalloc((void **)&d_pObjY, nPixYMap * sizeof(double));
  369. cudaMalloc((void **)&d_pObjZ, nSlices * sizeof(double));
  370. cudaCheckErrors("cudaMalloc Coordinates");
  371. // Pointer for mapped coordinates
  372. double* d_pDetmY;
  373. double* d_pDetmX;
  374. // Allocate global memory on the device for mapped coordinates
  375. cudaMalloc((void **)&d_pDetmY, nDetYMap * sizeof(double));
  376. cudaMalloc((void **)&d_pDetmX, nDetYMap * nDetXMap * sizeof(double));
  377. cudaCheckErrors("cudaMalloc Map-Coordinates");
  378. // Pointer for rotated detector coords
  379. double* d_pRdetY;
  380. double* d_pRdetZ;
  381. // Allocate global memory on the device for for rotated detector coords
  382. cudaMalloc((void **)&d_pRdetY, nDetYMap * sizeof(double));
  383. cudaMalloc((void **)&d_pRdetZ, nDetYMap * sizeof(double));
  384. cudaCheckErrors("cudaMalloc Rot-Coordinates");
  385. // Generate detector and object boudaries
  386. threadsPerBlock.x = maxThreadsPerBlock;
  387. blockSize.x = (nDetX / maxThreadsPerBlock) + 1;
  388. map_boudaries_kernel << <blockSize, threadsPerBlock >> > (d_pDetX, nDetXMap, (double)nDetX, -du, 0.0);
  389. blockSize.x = (nDetY / maxThreadsPerBlock) + 1;
  390. map_boudaries_kernel << <blockSize, threadsPerBlock >> > (d_pDetY, nDetYMap, nDetY / 2.0, dv, 0.0);
  391. blockSize.x = (nPixX / maxThreadsPerBlock) + 1;
  392. map_boudaries_kernel << <blockSize, threadsPerBlock >> > (d_pObjX, nPixXMap, (double)nPixX, -dx, 0.0);
  393. blockSize.x = (nPixY / maxThreadsPerBlock) + 1;
  394. map_boudaries_kernel << <blockSize, threadsPerBlock >> > (d_pObjY, nPixYMap, nPixY / 2.0, dy, 0.0);
  395. blockSize.x = (nSlices / maxThreadsPerBlock) + 1;
  396. map_boudaries_kernel << <blockSize, threadsPerBlock >> > (d_pObjZ, nSlices, 0.0, dz, DAG + (dz / 2.0));
  397. //mexPrintf("Map Boundaries -- Threads:%d Blocks:%d \n", threads, blocks);
  398. // Initiate variables value with 0
  399. cudaMemset(d_pDetZ, 0, nDetYMap * sizeof(double));
  400. cudaMemset(d_pProj, 0, nDetX * nDetY * nProj * sizeof(double));
  401. cudaCheckErrors("cudaMemset Zeros");
  402. //cudaDeviceSynchronize();
  403. // X - ray tube initial position
  404. double tubeX = 0;
  405. double tubeY = 0;
  406. double tubeZ = DSD;
  407. // Iso - center position
  408. double isoY = 0;
  409. double isoZ = DDR;
  410. // Integration of 2D slices over the whole volume
  411. // (S.1.Integration. - Liu et al(2017))
  412. // Naive integration o the X coord
  413. threadsPerBlock.x = 10;
  414. threadsPerBlock.y = 10;
  415. threadsPerBlock.z = 10;
  416. blockSize.x = (unsigned int) ceil(((double)nPixYMap / (threadsPerBlock.x-1)));
  417. blockSize.y = 1;
  418. blockSize.z = (unsigned int) ceil((double)nSlices / threadsPerBlock.z);
  419. //mexPrintf("Divisao: %f \n", (double)nPixY / (threadsPerBlock.x - 1));
  420. //mexPrintf("Divisao ceil: %d \n", (unsigned int)ceil((double)nPixY / (threadsPerBlock.x - 1)));
  421. for (int k = 0; k < ceil((double)nPixXMap / (threadsPerBlock.x-1)); k++) {
  422. img_integration_kernel << <blockSize, threadsPerBlock >> > (d_pVolume, nPixXMap, nPixYMap, integrateXcoord, 0, k * 9, nSlices);
  423. //mexPrintf("integration kernel -- threadsX:%d blocksX:%d threadsY:%d blocksY:%d \n", threadsPerBlock.x, blockSize.x, threadsPerBlock.y, blockSize.y);
  424. }
  425. // Naive integration o the Y coord
  426. threadsPerBlock.x = 10;
  427. threadsPerBlock.y = 10;
  428. threadsPerBlock.z = 10;
  429. blockSize.x = 1;
  430. blockSize.y = (unsigned int) ceil((double)nPixXMap / (threadsPerBlock.y-1));
  431. blockSize.z = (unsigned int) ceil((double)nSlices / threadsPerBlock.z);
  432. for (int k = 0; k < ceil((double)nPixYMap / (threadsPerBlock.y-1)); k++) {
  433. img_integration_kernel << <blockSize, threadsPerBlock >> > (d_pVolume, nPixXMap, nPixYMap, integrateYcoord, k * 9, 0, nSlices);
  434. //mexPrintf("integration kernel -- threadsX:%d blocksX:%d threadsY:%d blocksY:%d \n", threadsPerBlock.x, blockSize.x, threadsPerBlock.y, blockSize.y);
  435. }
  436. // Test if we will loop over all projs or not
  437. unsigned int projIni, projEnd;
  438. if(idXProj == -1){
  439. projIni = 0;
  440. projEnd = nProj;
  441. }
  442. else{
  443. projIni = (unsigned int) idXProj;
  444. projEnd = (unsigned int) idXProj + 1;
  445. }
  446. // For each projection
  447. for (unsigned int p = projIni; p < projEnd; p++) {
  448. // Get specif tube angle for the projection
  449. double theta = h_pTubeAngle[p] * M_PI / 180.0;
  450. // Get specif detector angle for the projection
  451. double phi = h_pDetAngle[p] * M_PI / 180.0;
  452. //mexPrintf("Tube angle:%f Det angle:%f\n", theta, phi);
  453. // Tube rotation
  454. double rtubeY = ((tubeY - isoY)*cos(theta) - (tubeZ - isoZ)*sin(theta)) + isoY;
  455. double rtubeZ = ((tubeY - isoY)*sin(theta) + (tubeZ - isoZ)*cos(theta)) + isoZ;
  456. //mexPrintf("R tube Y:%f R tube Z:%f\n", rtubeY, rtubeZ);
  457. // Detector rotation
  458. threadsPerBlock.x = maxThreadsPerBlock;
  459. threadsPerBlock.y = 1;
  460. threadsPerBlock.z = 1;
  461. blockSize.x = (nDetYMap / maxThreadsPerBlock) + 1;
  462. blockSize.y = 1;
  463. blockSize.z = 1;
  464. rot_detector_kernel << <blockSize, threadsPerBlock >> > (d_pRdetY, d_pRdetZ, d_pDetY, d_pDetZ, isoY, isoZ, phi, nDetYMap);
  465. //cudaDeviceSynchronize();
  466. //mexPrintf("Detector rotation -- Threads:%d Blocks:%d \n", threads, blocks);
  467. // For each slice
  468. for (unsigned int nz = 0; nz < nSlices; nz++) {
  469. /*
  470. Map detector onto XY plane(Inside proj loop in case detector rotates)
  471. *** Note: Matlab has linear indexing as a column of elements, i.e, the elements are actually stored in memory as queued columns.
  472. So threads in X are launched in the column direction (DBT Y coord), and threads in Y are launched in the row direction (DBT X coord).
  473. */
  474. threadsPerBlock.x = 32;
  475. threadsPerBlock.y = 32;
  476. threadsPerBlock.z = 1;
  477. blockSize.x = (nDetYMap / threadsPerBlock.x) + 1;
  478. blockSize.y = (nDetXMap / threadsPerBlock.y) + 1;
  479. blockSize.z = 1;
  480. mapDet2Slice_kernel << <blockSize, threadsPerBlock >> > (d_pDetmX, d_pDetmY, tubeX, rtubeY, rtubeZ, d_pDetX, d_pRdetY, d_pRdetZ, d_pObjZ, nDetXMap, nDetYMap, nz);
  481. //cudaDeviceSynchronize();
  482. //mexPrintf("Map detector onto XY plane -- ThreadsX:%d ThreadsY:%d BlocksX:%d BlocksY:%d\n", threadsPerBlock.x, threadsPerBlock.y, blockSizeX, blockSizeY);
  483. /*
  484. Creates texture memory
  485. */
  486. /*
  487. // Allocate CUDA array in device memory
  488. cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 0, 0, 0,cudaChannelFormatKindFloat);
  489. cudaArray* d_imageArray;
  490. cudaMallocArray(&d_imageArray, &channelDesc, nPixXMap, nPixYMap);
  491. // Copy to texture memory
  492. checkCudaErrors(cudaMemcpyToArray(d_imageArray, 0, 0, d_pSlice, nPixXMap* nPixYMap * sizeof(double), cudaMemcpyDeviceToDevice));
  493. // Specify texture
  494. struct cudaResourceDesc resDesc;
  495. memset(&resDesc, 0, sizeof(resDesc));
  496. resDesc.resType = cudaResourceTypeArray;
  497. resDesc.res.array.array = d_imageArray;
  498. // Specify texture object parameters
  499. struct cudaTextureDesc texDesc;
  500. memset(&texDesc, 0, sizeof(texDesc));
  501. texDesc.addressMode[0] = cudaAddressModeBorder;
  502. texDesc.addressMode[1] = cudaAddressModeBorder;
  503. texDesc.filterMode = cudaFilterModeLinear;
  504. texDesc.readMode = cudaReadModeElementType;
  505. texDesc.normalizedCoords = true;
  506. // Create texture object
  507. cudaTextureObject_t texObj = 0;
  508. cudaCreateTextureObject(&texObj, &resDesc, &texDesc, NULL);
  509. // Allocate result of transformation in device memory
  510. double* d_projI;
  511. cudaMalloc(&d_projI, nDetXMap*nDetYMap * sizeof(double));
  512. // Invoke kernel
  513. blockSize.x = (nDetYMap / threadsPerBlock.x) + 1;
  514. blockSize.y = (nDetXMap / threadsPerBlock.y) + 1;
  515. const double2 normFactor = make_double2(nPixX*-dx,
  516. nPixY*dy);
  517. const double2 offset = make_double2(1,
  518. 0.5);
  519. bilinear_interp_kernel << <blockSize, threadsPerBlock >> > (d_projI, texObj, d_pDetmX, d_pDetmY, nDetXMap, nDetYMap, normFactor, offset);
  520. mexPrintf("Texture interpolation -- Threads:%d Blocks:%d \n", threads, blocks);
  521. Destroy texture object
  522. cudaDestroyTextureObject(texObj);
  523. */
  524. /*
  525. S.2. Interpolation - Liu et al (2017)
  526. */
  527. blockSize.x = (nDetYMap / threadsPerBlock.x) + 1;
  528. blockSize.y = (nDetXMap / threadsPerBlock.y) + 1;
  529. bilinear_interpolation_kernel_GPU << <blockSize, threadsPerBlock >> > (d_projI, d_pVolume, d_pDetmX, d_pDetmY, nDetXMap, nDetYMap, nPixXMap, nPixYMap, nPixX, nPixY, dx, nz);
  530. /*
  531. S.3. Differentiation - Eq. 24 - Liu et al (2017)
  532. */
  533. blockSize.x = (nDetY / threadsPerBlock.x) + 1;
  534. blockSize.y = (nDetX / threadsPerBlock.y) + 1;
  535. differentiation_kernel << <blockSize, threadsPerBlock >> > (d_pProj, d_projI, d_pDetmX, d_pDetmY, tubeX, rtubeY, rtubeZ, d_pDetX, d_pRdetY, d_pRdetZ, nDetX, nDetY, nDetXMap, nDetYMap, du, dv, dx, dy, dz, p);
  536. } // Loop end slices
  537. } // Loop end Projections
  538. //nDetYMap
  539. //double* const h_pPixmX = (double*)mxMalloc(nDetXMap*nDetYMap * sizeof(double));
  540. //double* const h_pPixmY = (double*)mxMalloc(nDetYMap * sizeof(double));
  541. //cudaMemcpy((void *)h_pPixmX, (void *)d_pDetmX, nDetXMap*nDetYMap * sizeof(double), cudaMemcpyDeviceToHost);
  542. //cudaMemcpy((void *)h_pPixmY, (void *)d_pDetmY, nDetYMap * sizeof(double), cudaMemcpyDeviceToHost);
  543. //for (int u = 0; u < nDetXMap; u++)
  544. // for (int v = 0; v < nDetYMap; v++) {
  545. // h_pProj[(0 * nDetYMap*nDetXMap) + (u*nDetYMap) + v] = h_pPixmX[(u*nDetYMap) + v];
  546. // h_pProj[(1 * nDetYMap*nDetXMap) + (u*nDetYMap) + v] = h_pPixmY[v];
  547. // }
  548. // d_projI
  549. //double* const h_pPixmX = (double*)mxMalloc(nDetXMap*nDetYMap * sizeof(double));
  550. //cudaMemcpy((void *)h_pPixmX, (void *)d_projI, nDetXMap*nDetYMap * sizeof(double), cudaMemcpyDeviceToHost);
  551. //for (int u = 0; u < nDetXMap; u++)
  552. // for (int v = 0; v < nDetYMap; v++) {
  553. // h_pProj[(0 * nDetYMap*nDetXMap) + (u*nDetYMap) + v] = h_pPixmX[(u*nDetYMap) + v];
  554. // }
  555. // d_pVolume
  556. //cudaMemcpy((void *)h_pProj, (void *)d_pVolume, nSlices*nPixXMap*nPixYMap * sizeof(double), cudaMemcpyDeviceToHost);
  557. // d_pProj
  558. cudaMemcpy((void *)h_pProj, (void *)d_pProj, nProj* nDetX * nDetY * sizeof(double), cudaMemcpyDeviceToHost);
  559. cudaCheckErrors("cudaMemcpy Final");
  560. cudaFree(d_pProj);
  561. cudaFree(d_projI);
  562. cudaFree(d_pVolume);
  563. cudaFree(d_pTubeAngle);
  564. cudaFree(d_pDetAngle);
  565. cudaFree(d_pDetX);
  566. cudaFree(d_pDetY);
  567. cudaFree(d_pDetZ);
  568. cudaFree(d_pObjX);
  569. cudaFree(d_pObjY);
  570. cudaFree(d_pObjZ);
  571. cudaFree(d_pDetmY);
  572. cudaFree(d_pDetmX);
  573. cudaFree(d_pRdetY);
  574. cudaFree(d_pRdetZ);
  575. cudaCheckErrors("cudaFree Final");
  576. sdkStopTimer(&timer);
  577. mexPrintf("Processing time: %f (ms)\n", sdkGetTimerValue(&timer));
  578. sdkDeleteTimer(&timer);
  579. cudaDeviceReset();
  580. return;
  581. }