updatephi_kernel.cu
上传用户:liuping58
上传日期:2022-06-05
资源大小:105k
文件大小:3k
源码类别:

3D图形编程

开发平台:

Visual C++

  1. #define ALPHA  0.003
  2. #define DT  0.2
  3. #define max(x,y)    ((x>y) ? x : y )
  4. #define min(x,y)    ((x<y) ? x : y )
  5. __global__ void updatephi( float *d_phi, float *d_phi1, float *d_D, int imageW, int imageH)
  6. {
  7. int c= blockIdx.x * blockDim.x + threadIdx.x;
  8. int r= blockIdx.y * blockDim.y + threadIdx.y;
  9. int ind= r*imageW+c;
  10. if(ind<imageW*imageH){
  11. float dx,dxplus,dxminus,dxplusy,dxminusy;
  12. float dy, dyplus,dyminus,dyplusx,dyminusx;
  13. float gradphimax, gradphimin, nplusx, nplusy, nminusx, nminusy, curvature;
  14. float F, gradphi;
  15. if(c==0||c==imageW-1){dx=0;} else {dx=(d_phi1[ind+1]-d_phi1[ind-1])/2;}
  16. if(c==imageW-1){dxplus=0;} else {dxplus=(d_phi1[ind+1]-d_phi1[ind]);}
  17. if(c==0){dxminus=0;} else {dxminus=(d_phi1[ind]-d_phi1[ind-1]);}
  18. if(r==0||c==0||c==imageW-1){dxplusy=0;} else {dxplusy=(d_phi1[ind-imageW+1]-d_phi1[ind-imageW-1])/2;}
  19. if(r==imageH-1||c==0||c==imageW-1){dxminusy=0;} else {dxminusy=(d_phi1[ind+imageW+1]-d_phi1[ind+imageW-1])/2;}
  20. if(r==0||r==imageH-1){dy=0;} else {dy=(d_phi1[ind-imageW]-d_phi1[ind+imageW])/2;}
  21. if(r==0){dyplus=0;} else {dyplus=(d_phi1[ind-imageW]-d_phi1[ind]);}
  22. if(r==imageH-1){dyminus=0;} else {dyminus=(d_phi1[ind]-d_phi1[ind+imageW]);}
  23. if(r==0||c==imageW-1||r==imageH-1){dyplusx=0;} else {dyplusx=(d_phi1[ind-imageW+1]-d_phi1[ind+imageW+1])/2;}
  24. if(r==0||c==0||r==imageH-1){dyminusx=0;} else {dyminusx=(d_phi1[ind-imageW-1]-d_phi1[ind+imageW-1])/2;}
  25. 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)))
  26.    +(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))));
  27. 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)))
  28.    +(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))));
  29. nplusx= dxplus / sqrt(1.192092896e-07F + (dxplus*dxplus) + ((dyplusx + dy)*(dyplusx + dy)*0.25) );
  30. nplusy= dyplus / sqrt(1.192092896e-07F + (dyplus*dyplus) + ((dxplusy + dx)*(dxplusy + dx)*0.25) );
  31. nminusx= dxminus / sqrt(1.192092896e-07F + (dxminus*dxminus) + ((dyminusx + dy)*(dyminusx + dy)*0.25) );
  32. nminusy= dyminus / sqrt(1.192092896e-07F + (dyminus*dyminus) + ((dxminusy + dx)*(dxminusy + dx)*0.25) );
  33. curvature= ((nplusx-nminusx)+(nplusy-nminusy))/2;
  34. F = (-ALPHA * d_D[ind]) + ((1-ALPHA) * curvature);
  35. if(F>0) {gradphi=gradphimax;} else {gradphi=gradphimin;}
  36. d_phi[ind]=d_phi1[ind] + (DT * F * gradphi);
  37. }
  38. }