开通博客赚积分
发布资源赚积分
分类

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 3D图形编程 > 查看源码
updatephi_kernel.cu
资源名称：CUDAseg.zip　[点击查看]
上传用户：liuping58
上传日期：2022-06-05
资源大小：105k
文件大小：4k
源码类别：

3D图形编程

开发平台：

Visual C++

updatephi_kernel.cu：源码内容
  2. #define ALPHA  0.003
  3. #define DT  0.2
  5. #define max(x,y)    ((x>y) ? x : y )
  6. #define min(x,y)    ((x<y) ? x : y )
  8. #define INDEX(i,j,j_off)  (i +__mul24(j,j_off))
  10. #define BLOCKDIM_X  16
  11. #define BLOCKDIM_Y  32
  13. __global__ void updatephi( float *d_phi, float *d_phi1, float *d_D, int imageW, int imageH, int pitch)
  15. {
  16. int c= blockIdx.x * blockDim.x + threadIdx.x;
  17. int r= blockIdx.y * blockDim.y + threadIdx.y;
  18. int ind= r*pitch+c;
  20. int   indg, indg_h, indg0;
  21. int   i, j, k, ind2, ind_h, halo, active;
  23. __shared__ float s_data[(BLOCKDIM_X+2)*(BLOCKDIM_Y+2)];
  25. k    =  threadIdx.y*BLOCKDIM_X + threadIdx.x;
  26. halo = k < 2*(BLOCKDIM_X+BLOCKDIM_Y+2);
  28. if (halo) {
  29. if (threadIdx.y<2) {               // y-halos (coalesced)
  30. i = threadIdx.x;
  31. j = threadIdx.y*(BLOCKDIM_Y+1) - 1;
  32. }
  33. else {                             // x-halos (not coalesced)
  34. i = (k%2)*(BLOCKDIM_X+1) - 1;
  35. j =  k/2 - BLOCKDIM_X - 1;
  36. }
  38. ind_h  = INDEX(i+1,j+1,BLOCKDIM_X+2);
  40. i      = INDEX(i,blockIdx.x,BLOCKDIM_X);   // global indices
  41. j      = INDEX(j,blockIdx.y,BLOCKDIM_Y);
  42. indg_h = INDEX(i,j,pitch);
  44. halo   =  (i>=0) && (i<imageW) && (j>=0) && (j<imageH);
  45. }
  47. //
  48. // then set up indices for main block
  49. //
  51. i    = threadIdx.x;
  52. j    = threadIdx.y;
  53. ind2  = INDEX(i+1,j+1,BLOCKDIM_X+2) ;
  55. i    = INDEX(i,blockIdx.x,BLOCKDIM_X);     // global indices
  56. j    = INDEX(j,blockIdx.y,BLOCKDIM_Y);
  57. indg = INDEX(i,j,pitch);
  59. active = (i<imageW) && (j<imageH);
  61. //
  62. // read initial plane of u1 array
  63. //
  65. if (active) s_data[ind2] = d_phi1[indg];
  66. if (halo) s_data[ind_h] = d_phi1[indg_h];
  68. __syncthreads();
  71. if(active){
  73. if(r<imageW&&c<imageH){
  75. float dx,dxplus,dxminus,dxplusy,dxminusy;
  77. float dy, dyplus,dyminus,dyplusx,dyminusx;
  79. float gradphimax, gradphimin, nplusx, nplusy, nminusx, nminusy, curvature;
  80. float F, gradphi;
  82. if(c==0||c==imageW-1){dx=0;} else {dx=(s_data[ind2+1]-s_data[ind2-1])/2;}
  83. if(c==imageW-1){dxplus=0;} else {dxplus=(s_data[ind2+1]-s_data[ind2]);}
  84. if(c==0){dxminus=0;} else {dxminus=(s_data[ind2]-s_data[ind2-1]);}
  85. if(r==0||c==0||c==imageW-1){dxplusy=0;} else {dxplusy=(s_data[ind2-(BLOCKDIM_X+2)+1]-s_data[ind2-(BLOCKDIM_X+2)-1])/2;}
  86. if(r==imageH-1||c==0||c==imageW-1){dxminusy=0;} else {dxminusy=(s_data[ind2+(BLOCKDIM_X+2)+1]-s_data[ind2+(BLOCKDIM_X+2)-1])/2;}
  88. if(r==0||r==imageH-1){dy=0;} else {dy=(s_data[ind2-(BLOCKDIM_X+2)]-s_data[ind2+(BLOCKDIM_X+2)])/2;}
  89. if(r==0){dyplus=0;} else {dyplus=(s_data[ind2-(BLOCKDIM_X+2)]-s_data[ind2]);}
  90. if(r==imageH-1){dyminus=0;} else {dyminus=(s_data[ind2]-s_data[ind2+(BLOCKDIM_X+2)]);}
  91. if(r==0||c==imageW-1||r==imageH-1){dyplusx=0;} else {dyplusx=(s_data[ind2-(BLOCKDIM_X+2)+1]-s_data[ind2+(BLOCKDIM_X+2)+1])/2;}
  92. if(r==0||c==0||r==imageH-1){dyminusx=0;} else {dyminusx=(s_data[ind2-(BLOCKDIM_X+2)-1]-s_data[ind2+(BLOCKDIM_X+2)-1])/2;}
  94. gradphimax=sqrt((sqrt(max(dxplus,0)*max(dxplus,0)+max(-dxminus,0)*max(-dxminus,0)))*(sqrt(max(dxplus,0)*max(dxplus,0)+max(-dxminus,0)*max(-dxminus,0)))
  95. +(sqrt(max(dyplus,0)*max(dyplus,0)+max(-dyminus,0)*max(-dyminus,0)))*(sqrt(max(dyplus,0)*max(dyplus,0)+max(-dyminus,0)*max(-dyminus,0))));
  97. gradphimin=sqrt((sqrt(min(dxplus,0)*min(dxplus,0)+min(-dxminus,0)*min(-dxminus,0)))*(sqrt(min(dxplus,0)*min(dxplus,0)+min(-dxminus,0)*min(-dxminus,0)))
  98. +(sqrt(min(dyplus,0)*min(dyplus,0)+min(-dyminus,0)*min(-dyminus,0)))*(sqrt(min(dyplus,0)*min(dyplus,0)+min(-dyminus,0)*min(-dyminus,0))));
  100. nplusx= dxplus / sqrt(1.192092896e-07F + (dxplus*dxplus) + ((dyplusx + dy)*(dyplusx + dy)*0.25) );
  101. nplusy= dyplus / sqrt(1.192092896e-07F + (dyplus*dyplus) + ((dxplusy + dx)*(dxplusy + dx)*0.25) );
  102. nminusx= dxminus / sqrt(1.192092896e-07F + (dxminus*dxminus) + ((dyminusx + dy)*(dyminusx + dy)*0.25) );
  103. nminusy= dyminus / sqrt(1.192092896e-07F + (dyminus*dyminus) + ((dxminusy + dx)*(dxminusy + dx)*0.25) );
  104. curvature= ((nplusx-nminusx)+(nplusy-nminusy))/2;
  106. F = (-ALPHA * d_D[indg]) + ((1-ALPHA) * curvature);
  107. if(F>0) {gradphi=gradphimax;} else {gradphi=gradphimin;}
  108. d_phi[indg]=s_data[ind2] + (DT * F * gradphi);
  110. }
  112. __syncthreads();
  117. }
  118. }