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

源码开发语言/平台
当前位置：首页 > 源码/资料 > 图形图像 > 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 <cutil.h>
  6. #include <math.h>
  9. #define max(x,y)    ((x>y) ? x : y )
  10. #define min(x,y)    ((x<y) ? x : y )
  12. #define DT 0.15
  14. char *volumeFilename, *maskFilename;
  15. int ITERATIONS, THRESHOLD, EPSILON;
  16. float alpha;
  18. float *phi, *D, *contour;
  19. size_t size;
  20. unsigned char *input,*output;
  21. int imageW, imageH, imageD, N;
  23. int its=0;
  24. unsigned int Timer = 0;
  25. unsigned int IterationTimer = 0;
  27. int i,j,k;
  31. unsigned char* loadRawUchar(char *filename, size_t size){
  33. FILE *fp = fopen(filename, "rb");
  34.     if (!fp) {
  35.         fprintf(stderr, "Error opening file '%s'n", filename);
  36.         return 0;
  37.     }
  39. unsigned char *data = (unsigned char *) malloc(size);
  40. size_t read = fread(data, 1, size, fp);
  41. fclose(fp);
  43.     printf("Read '%s', %d bytesn", filename, read);
  45.     return data;
  46. }
  48. float *loadMask(char *filename, size_t size){
  50. FILE *fp = fopen(filename, "rb");
  51.     if (!fp) {
  52.         fprintf(stderr, "Error opening file '%s'n", filename);
  53.         return 0;
  54.     }
  56. float *data = (float *) malloc(size*sizeof(float));
  57. size_t read = fread(data, 1, size, fp);
  58. fclose(fp);
  60.     printf("Read '%s', %d bytesn", filename, read);
  62.     return data;
  63. }
  65. void *writeoutput(unsigned char *data, size_t size){
  67. char *outputFilename= "output.raw";
  68. FILE *fp = fopen(outputFilename, "w");
  69.     if (!fp) {
  70.         fprintf(stderr, "Error opening file '%s'n", outputFilename);
  71.         return 0;
  72.     }
  74. size_t write = fwrite(data, 1, size, fp);
  75. fclose(fp);
  77.     printf("Write '%s', %d bytesn", outputFilename, write);
  79.     return 0;
  80. }
  82. void update_phi(){
  84. float dx,dy,dz;
  85. float dxplus, dyplus, dzplus, dxminus, dyminus, dzminus;
  86. float dxplusy, dxminusy, dxplusz, dxminusz, dyplusx, dyminusx, dyplusz, dyminusz, dzplusx, dzminusx, dzplusy, dzminusy;
  88. float gradphimax, gradphimin, nplusx, nplusy, nplusz, nminusx, nminusy, nminusz, curvature;
  89. float F, gradphi;
  91. float *phi1;
  93. int ind, jOFF, kOFF;
  95. jOFF=imageW;
  96. kOFF=imageW*imageH;
  98. phi1=(float *)malloc(imageW*imageH*imageD*sizeof(float));
  99. for(i=0;i<N;i++){
  100. phi1[i]=phi[i];
  101. }
  103. for(i=0;i<imageW;i++){
  104. for(j=0;j<imageH;j++){
  105. for(k=0;k<imageD;k++){
  107. ind=i+j*imageW+k*imageW*imageH;
  109. if(i==0||i==imageW-1){dx=0;} else {dx=(phi1[ind+1]-phi1[ind-1])/2;}
  110. if(j==0||j==imageH-1){dy=0;} else {dy=(phi1[ind-imageW]-phi1[ind+imageW])/2;}
  111. if(k==0||k==imageD-1){dz=0;} else {dz=(phi1[ind+kOFF]-phi1[ind-kOFF])/2;}
  113. if(i==imageW-1){dxplus=0;} else {dxplus=(phi1[ind+1]-phi1[ind]);}
  114. if(j==0){dyplus=0;} else {dyplus=(phi1[ind-imageW]-phi1[ind]);}
  115. if(k==imageD-1){dzplus=0;} else {dzplus=(phi1[ind+kOFF]-phi1[ind]);}
  116. if(i==0){dxminus=0;} else {dxminus=(phi1[ind]-phi1[ind-1]);}
  117. if(j==imageH-1){dyminus=0;} else {dyminus=(phi1[ind]-phi1[ind+imageW]);}
  118. if(k==0){dzminus=0;} else {dzminus=(phi1[ind]-phi1[ind-kOFF]);}
  120. if(i==0||i==imageW-1||j==0){dxplusy=0;} else {dxplusy=(phi1[ind-imageW+1]-phi1[ind-imageW-1])/2;}
  121. if(i==0||i==imageW-1||j==imageH-1){dxminusy=0;} else {dxminusy=(phi1[ind+imageW+1]-phi1[ind+imageW-1])/2;}
  122. if(i==0||i==imageW-1||k==imageD-1) {dxplusz=0;} else {dxplusz=(phi1[ind+kOFF+1]-phi1[ind+kOFF-1])/2;}
  123. if(i==0||i==imageW-1||k==0) {dxminusz=0;} else {dxminusz=(phi1[ind-kOFF+1]-phi1[ind-kOFF-1])/2;}
  124. if(j==0||j==imageH-1||i==imageW-1){dyplusx=0;} else {dyplusx=(phi1[ind-imageW+1]-phi1[ind+imageW+1])/2;}
  125. if(j==0||j==imageH-1||i==0){dyminusx=0;} else {dyminusx=(phi1[ind-imageW-1]-phi1[ind+imageW-1])/2;}
  126. if(j==0||j==imageH-1||k==imageD-1) {dyplusz=0;} else {dyplusz=(phi1[ind+kOFF-jOFF]-phi1[ind+kOFF+jOFF])/2;}
  127. if(j==0||j==imageH-1||k==0) {dyminusz=0;} else {dyminusz=(phi1[ind-kOFF-jOFF]-phi1[ind-kOFF+jOFF])/2;}
  128. if(k==0||k==imageD-1||i==imageW-1) {dzplusx=0;} else {dzplusx=(phi1[ind+1+kOFF]-phi1[ind+1-kOFF])/2;}
  129. if(k==0||k==imageD-1||i==0) {dzminusx=0;} else {dzminusx=(phi1[ind-1+kOFF]-phi1[ind-1-kOFF])/2;}
  130. if(k==0||k==imageD-1||j==0) {dzplusy=0;} else {dzplusy=(phi1[ind-jOFF+kOFF]-phi1[ind-jOFF-kOFF])/2;}
  131. if(k==0||k==imageD-1||j==imageH-1) {dzminusy=0;} else {dzminusy=(phi1[ind+jOFF+kOFF]-phi1[ind+jOFF-kOFF])/2;}
  134. 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)))
  135. +(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)))
  136. +(sqrt(max(dzplus,0)*max(dzplus,0)+max(-dzminus,0)*max(-dzminus,0)))*(sqrt(max(dzplus,0)*max(dzplus,0)+max(-dzminus,0)*max(-dzminus,0))));
  138. 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)))
  139. +(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)))
  140. +(sqrt(min(dzplus,0)*min(dzplus,0)+min(-dzminus,0)*min(-dzminus,0)))*(sqrt(min(dzplus,0)*min(dzplus,0)+min(-dzminus,0)*min(-dzminus,0))));
  142. nplusx= dxplus / sqrt(1.192092896e-07F + (dxplus*dxplus) + ((dyplusx + dy)*(dyplusx + dy)*0.25) + ((dzplusx + dz)*(dzplusx + dz)*0.25));
  143. nplusy= dyplus / sqrt(1.192092896e-07F + (dyplus*dyplus) + ((dxplusy + dx)*(dxplusy + dx)*0.25) + ((dzplusy + dz)*(dzplusy + dz)*0.25));
  144. nplusz= dzplus / sqrt(1.192092896e-07F + (dzplus*dzplus) + ((dxplusz + dz)*(dxplusz + dz)*0.25) + ((dyplusz + dy)*(dyplusz + dy)*0.25));
  146. nminusx= dxminus / sqrt(1.192092896e-07F + (dxminus*dxminus) + ((dyminusx + dy)*(dyminusx + dy)*0.25) + ((dzminusx + dz)*(dzminusx + dz)*0.25));
  147. nminusy= dyminus / sqrt(1.192092896e-07F + (dyminus*dyminus) + ((dxminusy + dx)*(dxminusy + dx)*0.25) + ((dzminusy + dz)*(dzminusy + dz)*0.25));
  148. nminusz= dzminus / sqrt(1.192092896e-07F + (dzminus*dzminus) + ((dxminusz + dz)*(dxminusz + dz)*0.25) + ((dyminusz + dy)*(dyminusz + dy)*0.25));
  150. curvature= ((nplusx-nminusx)+(nplusy-nminusy)+(nplusz-nminusz))/2;
  152. F = (-alpha * D[ind]) + ((1-alpha) * curvature);
  153. if(F>0) {gradphi=gradphimax;} else {gradphi=gradphimin;}
  154. phi[ind]=phi1[ind] + (DT * F * gradphi);
  156. }
  157. }
  158. }
  161. free(phi1);
  163. }
  167. int main(int argc, char** argv){
  169. if(argc<9){
  170. printf("Too few command line arguments specified. Example: Seg -volume=brain_181_217_181.raw -mask=phi.raw -xsize=181 -ysize=217 -zsize=181 -iterations=1000 -threshold=150 -epsilon=50 -alpha=0.01n");
  171. exit(0);
  172. }
  174. cutGetCmdLineArgumentstr( argc, (const char**) argv, "volume", &volumeFilename);
  175. cutGetCmdLineArgumentstr( argc, (const char**) argv, "mask", &maskFilename);
  176. cutGetCmdLineArgumenti( argc, (const char**) argv, "xsize", &imageW);
  177. cutGetCmdLineArgumenti( argc, (const char**) argv, "ysize", &imageH);
  178. cutGetCmdLineArgumenti( argc, (const char**) argv, "zsize", &imageD);
  179. cutGetCmdLineArgumenti( argc, (const char**) argv, "iterations", &ITERATIONS);
  180. cutGetCmdLineArgumenti( argc, (const char**) argv, "threshold", &THRESHOLD);
  181. cutGetCmdLineArgumenti( argc, (const char**) argv, "epsilon", &EPSILON);
  182. cutGetCmdLineArgumentf( argc, (const char**) argv, "alpha", &alpha);
  184. N=imageW*imageH*imageD;
  185. input = loadRawUchar(volumeFilename, N);
  186. phi = loadMask(maskFilename, N);
  188. if((D = (float *)malloc(imageW*imageH*imageD*sizeof(float)))==NULL)printf("ME_Dn");
  189. for(i=0;i<N;i++){
  190. D[i] = EPSILON - abs(input[i] - THRESHOLD);
  191. }
  193. // Set up CUDA Timer
  194. cutCreateTimer(&Timer);
  195. cutCreateTimer(&IterationTimer);
  197. cutStartTimer(Timer);
  199. for(its=0;its<ITERATIONS;its++){
  200. update_phi();
  201. if(its%10==0)printf("Iteration %3d Total Time: %3.2f ReInit Time: %3.2fn", its, 0.001*cutGetTimerValue(Timer), 0.001*cutGetTimerValue(IterationTimer));
  202. }
  204. if((output = (unsigned char *) malloc(size))==NULL)printf("ME_OUTPUTn");
  205. for(i=0;i<N;i++){
  206. if(phi[i]>0){output[i]=0;} else { output[i]=255; }
  207. }
  209. writeoutput(output, size);
  213. }