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

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

3D图形编程

开发平台：

Visual C++

main.cpp：源码内容
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <float.h>
  5. #include <GL/glew.h>
  6. #include <cuda_runtime.h>
  7. #include <cutil.h>
  8. #include <math.h>
  9. #include <GL/glut.h>
  11. #define IMAGE "bigbrain.bmp"
  12. #define MASK "bigmask.bmp"
  15. #define ITERATIONS   5000
  16. #define THRESHOLD  170
  17. #define EPSILON  35
  19. #define ALPHA  0.009
  20. #define DT  0.25
  22. #define RITS  50
  24. float *phi, *phi1, *D, *contour;
  25. uchar4 *h_Src, *h_Mask;
  26. int imageW, imageH, N;
  30. void LoadBMPFile(uchar4 **dst, int *width, int *height, const char *name);
  31. void sedt2d(int *_d,unsigned char *_bimg,int _h,int _w);
  34. int its=0;
  35. unsigned int Timer = 0;
  36. unsigned int ReInitTimer = 0;
  38. int r;
  39. int c;
  40. int i;
  42. void init_phi(){
  44. int *init;
  45. unsigned char *mask;
  46. const char *mask_path = MASK;
  47. if((init=(int *)malloc(imageW*imageH*sizeof(int)))==NULL)printf("ME_INITn");
  48. if((phi=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_PHIn");
  49. mask = (unsigned char *)malloc(imageW*imageH*sizeof(unsigned char));
  51. //printf("Init Maskn");
  52. LoadBMPFile(&h_Mask, &imageW, &imageH, mask_path);
  55. for(r=0;r<imageH;r++){
  56. for(c=0;c<imageW;c++){
  57. mask[r*imageW+c] = (h_Mask[r*imageW+c].x)/255;
  58. //printf("%3d ", mask[r*imageW+c]);
  59. }
  60. //printf("n");
  61. }
  65. sedt2d(init,mask,imageH,imageW);
  69. //printf("sdf of init maskn");
  70. for(r=0;r<imageH;r++){
  71. for(c=0;c<imageW;c++){
  72. phi[r*imageW+c]=(float)init[r*imageW+c];
  73. if(phi[r*imageW+c]>0){
  74. phi[r*imageW+c]=0.5*sqrt(abs(phi[r*imageW+c]));
  75. } else {
  76. phi[r*imageW+c]=-0.5*sqrt(abs(phi[r*imageW+c]));
  77. }
  78. //printf("%6.3f ", phi[r*imageW+c]);
  79. }
  80. //printf("n");
  81. }
  83. free(init);
  84. free(mask);
  85. }
  87. void reinit_phi(){
  89. int *intphi;
  90. unsigned char *reinit;
  91. if((intphi=(int *)malloc(imageW*imageH*sizeof(int)))==NULL)printf("ME_INITn");
  92. reinit=(unsigned char *)malloc(imageW*imageH*sizeof(unsigned char));//TODO check
  94. for(i=0;i<N;i++){
  95. if(phi[i]<0){
  96. phi[i]=1;
  97. } else {
  98. phi[i]=0;
  99. }
  100. reinit[i]=(int)phi[i];
  101. }
  104. sedt2d(intphi,reinit,imageH,imageW);
  106. for(r=0;r<imageH;r++){
  107. for(c=0;c<imageW;c++){
  108. phi[r*imageW+c]=(float)intphi[r*imageW+c];
  109. if(phi[r*imageW+c]>0){
  110. phi[r*imageW+c]=0.5*sqrt(abs(phi[r*imageW+c]));
  111. } else {
  112. phi[r*imageW+c]=-0.5*sqrt(abs(phi[r*imageW+c]));
  113. }
  114. //printf("%6.3f ", phi[r*imageW+c]);
  115. }
  116. //printf("n");
  117. }
  119. free(reinit);
  120. free(intphi);
  121. }
  124. void update_phi(){
  126. float dx;
  127. float dxplus;
  128. float dxminus;
  129. float dxplusy;
  130. float dxminusy;
  131. float maxdxplus;
  132. float maxminusdxminus;
  133. float mindxplus;
  134. float minminusdxminus;
  136. float dy;
  137. float dyplus;
  138. float dyminus;
  139. float dyplusx;
  140. float dyminusx;
  141. float maxdyplus;
  142. float maxminusdyminus;
  143. float mindyplus;
  144. float minminusdyminus;
  146. float gradphimax;
  147. float gradphimin;
  149. float nplusx;
  150. float nplusy;
  151. float nminusx;
  152. float nminusy;
  153. float curvature;
  155. float F;
  156. float gradphi;
  159. int ind;
  163. for(i=0;i<N;i++){
  164. phi1[i]=phi[i];
  165. }
  167. for(r=0;r<imageH;r++){
  168. for(c=0;c<imageW;c++){
  170. ind=r*imageW+c;
  172. if(c==0||c==imageW-1){dx=0;} else {dx=(phi1[ind+1]-phi1[ind-1])/2;}
  173. if(c==imageW-1){dxplus=0;} else {dxplus=(phi1[ind+1]-phi1[ind]);}
  174. if(c==0){dxminus=0;} else {dxminus=(phi1[ind]-phi1[ind-1]);}
  175. if(r==0||c==0||c==imageW-1){dxplusy=0;} else {dxplusy=(phi1[ind-imageW+1]-phi1[ind-imageW-1])/2;}
  176. if(r==imageH-1||c==0||c==imageW-1){dxminusy=0;} else {dxminusy=(phi1[ind+imageW+1]-phi1[ind+imageW-1])/2;}
  177. if(dxplus<0){maxdxplus=0;} else { maxdxplus= dxplus*dxplus; }
  178. if(-dxminus<0){maxminusdxminus=0;} else { maxminusdxminus= dxminus*dxminus; }
  179. if(dxplus>0){mindxplus=0;} else { mindxplus= dxplus*dxplus; }
  180. if(-dxminus>0){minminusdxminus=0;} else { minminusdxminus= dxminus*dxminus; }
  182. if(r==0||r==imageH-1){dy=0;} else {dy=(phi1[ind-imageW]-phi1[ind+imageW])/2;}
  183. if(r==0){dyplus=0;} else {dyplus=(phi1[ind-imageW]-phi1[ind]);}
  184. if(r==imageH-1){dyminus=0;} else {dyminus=(phi1[ind]-phi1[ind+imageW]);}
  185. if(r==0||c==imageW-1||r==imageH-1){dyplusx=0;} else {dyplusx=(phi1[ind-imageW+1]-phi1[ind+imageW+1])/2;}
  186. if(r==0||c==0||r==imageH-1){dyminusx=0;} else {dyminusx=(phi1[ind-imageW-1]-phi1[ind+imageW-1])/2;}
  187. if(dyplus<0){maxdyplus=0;} else { maxdyplus= dyplus*dyplus; }
  188. if(-dyminus<0){maxminusdyminus=0;} else { maxminusdyminus= dyminus*dyminus; }
  189. if(dyplus>0){mindyplus=0;} else { mindyplus= dyplus*dyplus; }
  190. if(-dyminus>0){minminusdyminus=0;} else { minminusdyminus= dyminus*dyminus; }
  192. gradphimax=sqrt((sqrt(maxdxplus+maxminusdxminus))*(sqrt(maxdxplus+maxminusdxminus))+(sqrt(maxdyplus+maxminusdyminus))*(sqrt(maxdyplus+maxminusdyminus)));
  193. gradphimin=sqrt((sqrt(mindxplus+minminusdxminus))*(sqrt(mindxplus+minminusdxminus))+(sqrt(mindyplus+minminusdyminus))*(sqrt(mindyplus+minminusdyminus)));
  194. nplusx= dxplus / sqrt(1.192092896e-07F + (dxplus*dxplus) + ((dyplusx + dy)*(dyplusx + dy)*0.25) );
  195. nplusy= dyplus / sqrt(1.192092896e-07F + (dyplus*dyplus) + ((dxplusy + dx)*(dxplusy + dx)*0.25) );
  196. nminusx= dxminus / sqrt(1.192092896e-07F + (dxminus*dxminus) + ((dyminusx + dy)*(dyminusx + dy)*0.25) );
  197. nminusy= dyminus / sqrt(1.192092896e-07F + (dyminus*dyminus) + ((dxminusy + dx)*(dxminusy + dx)*0.25) );
  198. curvature= ((nplusx-nminusx)+(nplusy-nminusy))/2;
  200. F = (-ALPHA * D[ind]) + ((1-ALPHA) * curvature);
  201. if(F>0) {gradphi=gradphimax;} else {gradphi=gradphimin;}
  202. phi[ind]=phi1[ind] + (DT * F * gradphi);
  203. }
  204. }
  208. }
  212. void disp(void){
  215. glClear(GL_COLOR_BUFFER_BIT);
  218. update_phi();
  221. its++;
  223. if(its<ITERATIONS){
  224. glutPostRedisplay();
  226. if(its%50==0){
  228. printf("Iteration %3d Total Time: %3.2f ReInit Time: %3.2fn", its, 0.001*cutGetTimerValue(Timer), 0.001*cutGetTimerValue(ReInitTimer));
  230. cutStartTimer(ReInitTimer); // ReInit Timer Start
  232. reinit_phi(); // ReInit
  234. glDrawPixels(imageW, imageH, GL_GREEN, GL_FLOAT, phi);
  235. glutSwapBuffers();
  236. cutStopTimer(ReInitTimer); // ReInit Timer Stop
  237. }
  239. } else {
  241. printf("Iteration %3d Total Time: %3.2f ReInit Time: %3.2fn", its, 0.001*cutGetTimerValue(Timer), 0.001*cutGetTimerValue(ReInitTimer));
  243. glDrawPixels(imageW, imageH, GL_GREEN, GL_FLOAT, phi);
  244. glutSwapBuffers();
  248. }
  250. }
  252. int main(int argc, char** argv){
  254. const char *image_path = IMAGE;
  256. //TODO : declare ALL variables here
  258. LoadBMPFile(&h_Src, &imageW, &imageH, image_path);
  259. D = (float *)malloc(imageW*imageH*sizeof(float));
  261. //printf("Input Imagen");
  262. for(r=0;r<imageH;r++){
  263. for(c=0;c<imageW;c++){
  264. D[r*imageW+c] = h_Src[r*imageW+c].x;
  265. /*printf("%3.0f ", D[r*imageW+c]);*/
  266. }
  267. //printf("n");
  268. }
  270. N = imageW*imageH;
  272. for(i=0;i<N;i++){
  273. D[i] = EPSILON - abs(D[i] - THRESHOLD);
  274. }
  276. //printf("Speed Functionn");
  277. //for(int r=0;r<imageH;r++){
  278. // for(int c=0;c<imageW;c++){
  279. // printf("%3.0f ", D[r*imageW+c]);
  280. // }
  281. // printf("n");
  282. //}
  284. // Set up CUDA Timer
  285. cutCreateTimer(&Timer);
  286. cutCreateTimer(&ReInitTimer);
  288. cutStartTimer(Timer);
  290. init_phi();
  291. if((contour=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("Contourn");
  292. if((phi1=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("GRADPHIn");
  293. //update_phi();
  295.   // GL initialisation
  296.   glutInit(&argc, argv);
  297.   glutInitDisplayMode(GLUT_ALPHA | GLUT_DOUBLE);
  298.   glutInitWindowSize(imageW,imageH);
  299.   glutInitWindowPosition(100,100);
  300.   glutCreateWindow("GL Level Set Evolution");
  301.   glClearColor(0.0,0.0,0.0,0.0);
  304.   glutDisplayFunc(disp);
  305.   glutMainLoop();
  307. //printf("phi+1n");
  308. //for(int r=0;r<imageH;r++){
  309. // for(int c=0;c<imageW;c++){
  310. // printf("%6.3f ", phi[r*imageW+c]);
  311. // }
  312. // printf("n");
  313. //}
  314. }
  317. //TODO Memory Malloc Free
  319. //TODO Timer
  321. //TODO Comment Code