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

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

3D图形编程

开发平台：

Visual C++

main.cpp：源码内容
  1. #include <stdio.h>
  2. #include <stdlib.h>
  3. #include <string.h>
  4. #include <float.h>
  5. #include <cuda_runtime.h>
  6. #include <cutil.h>
  7. #include <math.h>
  8. #include <GL/glut.h>
  10. #define IMAGE "liver.bmp"
  12. #define ITERATIONS   5000
  13. #define THRESHOLD  170
  14. #define EPSILON  35
  16. #define ALPHA  0.009
  17. #define DT  0.25
  19. #define RITS  50
  21. float *phi, *D, *contour;
  22. uchar4 *h_Src, *h_Mask;
  23. int imageW, imageH, N;
  24. unsigned char *mask;
  27. void LoadBMPFile(uchar4 **dst, int *width, int *height, const char *name);
  28. void sedt2d(int *_d,unsigned char *_bimg,int _h,int _w);
  31. int its=0;
  32. unsigned int Timer = 0;
  35. void init_phi(){
  37. const char *mask_path = "mask.bmp";
  39. //printf("Init Maskn");
  40. LoadBMPFile(&h_Mask, &imageW, &imageH, mask_path);
  41. mask = (unsigned char *)malloc(imageW*imageH*sizeof(unsigned char));
  43. for(int r=0;r<imageH;r++){
  44. for(int c=0;c<imageW;c++){
  45. mask[r*imageW+c] = (h_Mask[r*imageW+c].x)/255;
  46. //printf("%3d ", mask[r*imageW+c]);
  47. }
  48. //printf("n");
  49. }
  51. int *init;
  52. if((init=(int *)malloc(imageW*imageH*sizeof(int)))==NULL)printf("ME_INITn");
  53. if((phi=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_PHIn");
  54. sedt2d(init,mask,imageH,imageW);
  58. //printf("sdf of init maskn");
  59. for(int r=0;r<imageH;r++){
  60. for(int c=0;c<imageW;c++){
  61. phi[r*imageW+c]=(float)init[r*imageW+c];
  62. if(phi[r*imageW+c]>0){
  63. phi[r*imageW+c]=0.5*sqrt(abs(phi[r*imageW+c]));
  64. } else {
  65. phi[r*imageW+c]=-0.5*sqrt(abs(phi[r*imageW+c]));
  66. }
  67. //printf("%6.3f ", phi[r*imageW+c]);
  68. }
  69. //printf("n");
  70. }
  72. free(init);
  73. free(mask);
  74. }
  76. void reinit_phi(){
  78. unsigned char *reinit;
  79. reinit=(unsigned char *)malloc(imageW*imageH*sizeof(unsigned char));
  80. int *intphi;
  81. if((intphi=(int *)malloc(imageW*imageH*sizeof(int)))==NULL)printf("ME_INITn");
  83. for(int i=0;i<N;i++){
  84. if(phi[i]<0){
  85. phi[i]=1;
  86. } else {
  87. phi[i]=0;
  88. }
  89. reinit[i]=(int)phi[i];
  90. }
  93. sedt2d(intphi,reinit,imageH,imageW);
  95. printf("reinitn");
  96. for(int r=0;r<imageH;r++){
  97. for(int c=0;c<imageW;c++){
  98. phi[r*imageW+c]=(float)intphi[r*imageW+c];
  99. if(phi[r*imageW+c]>0){
  100. phi[r*imageW+c]=0.5*sqrt(abs(phi[r*imageW+c]));
  101. } else {
  102. phi[r*imageW+c]=-0.5*sqrt(abs(phi[r*imageW+c]));
  103. }
  104. //printf("%6.3f ", phi[r*imageW+c]);
  105. }
  106. //printf("n");
  107. }
  109. free(reinit);
  110. free(intphi);
  111. }
  114. void update_phi(){
  116. float *ptr2phi;
  118. float *dx, *ptr2dx;
  119. if((dx=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_DXn");
  120. ptr2dx=dx;
  121. ptr2phi=phi;
  122. for(int r=0;r<imageH;r++){
  123. *ptr2dx=0;
  124. for(int c=0;c<imageW-2;c++){
  125. ptr2phi++;
  126. ptr2dx++;
  127. *ptr2dx = (*(ptr2phi+1) - *(ptr2phi-1)) /2 ;
  128. }
  129. ptr2dx++;
  130. *ptr2dx=0;
  131. ptr2dx++;
  132. ptr2phi++;
  133. if(r!=(imageH-1))ptr2phi++;
  134. }
  136. float *dxplus, *ptr2dxplus;
  137. if((dxplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_DX+n");
  138. for(int i=0;i<N;i++)dxplus[i]=phi[i];
  139. ptr2dxplus=dxplus;
  140. for(int r=0;r<imageH;r++){
  141. for(int c=0;c<imageW-1;c++){
  142. *ptr2dxplus = *(ptr2dxplus+1) - *ptr2dxplus;
  143. ptr2dxplus++;
  144. }
  145. *ptr2dxplus=0;
  146. ptr2dxplus++;
  147. }
  149. float *dxminus, *ptr2dxminus;
  150. if((dxminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_DX-n");
  151. for(int i=0;i<N;i++)dxminus[i]=phi[i];
  152. ptr2dxminus=dxminus+N-1;
  153. for(int r=0;r<imageH;r++){
  154. for(int c=0;c<imageW-1;c++){
  155. *ptr2dxminus = *ptr2dxminus - *(ptr2dxminus-1);
  156. ptr2dxminus--;
  157. }
  158. *ptr2dxminus=0;
  159. ptr2dxminus--;
  161. }
  163. float *dxplusy, *ptr2dxplusy;
  164. if((dxplusy=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dxplusyn");
  165. ptr2dxplusy=dxplusy;
  166. ptr2phi=phi;
  167. for(int c=0;c<imageW;c++){
  168. *ptr2dxplusy = 0;
  169. ptr2dxplusy++;
  170. }
  171. for(int r=0;r<imageH-1;r++){
  172. *ptr2dxplusy=0;
  173. for(int c=0;c<imageW-2;c++){
  174. ptr2phi++;
  175. ptr2dxplusy++;
  176. *ptr2dxplusy = (*(ptr2phi+1) - *(ptr2phi-1)) /2 ;
  177. }
  178. ptr2dxplusy++;
  179. *ptr2dxplusy=0;
  180. ptr2dxplusy++;
  181. ptr2phi++;
  182. if(r!=(imageH-2))ptr2phi++;
  183. }
  185. float *dxminusy, *ptr2dxminusy;
  186. if((dxminusy=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dxminusyn");
  187. ptr2dxminusy=dxminusy;
  188. ptr2phi=phi+imageW;
  189. for(int r=0;r<imageH-1;r++){
  190. *ptr2dxminusy=0;
  191. for(int c=0;c<imageW-2;c++){
  192. ptr2phi++;
  193. ptr2dxminusy++;
  194. *ptr2dxminusy = (*(ptr2phi+1) - *(ptr2phi-1) );
  195. }
  196. ptr2dxminusy++;
  197. *ptr2dxminusy=0;
  198. ptr2dxminusy++;
  199. ptr2phi++;
  200. if(r!=(imageH-2))ptr2phi++;
  201. }
  202. for(int c=0;c<imageW;c++){
  203. *ptr2dxminusy = 0;
  204. ptr2dxminusy++;
  205. }
  207. float *maxdxplus, *ptr2maxdxplus;
  208. if((maxdxplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_MDX+n");
  209. for(int i=0;i<N;i++)maxdxplus[i]=dxplus[i];
  210. ptr2maxdxplus = maxdxplus;
  211. for(int i=0;i<N;i++){
  212. if (*ptr2maxdxplus < 0) {
  213. *ptr2maxdxplus = 0;
  214. } else {
  215. *ptr2maxdxplus *= *ptr2maxdxplus;
  216. }
  217. ptr2maxdxplus++;
  218. }
  220. float *maxminusdxminus, *ptr2maxminusdxminus;
  221. if((maxminusdxminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_-MDX-n");
  222. for(int i=0;i<N;i++)maxminusdxminus[i]=-dxminus[i];
  223. ptr2maxminusdxminus = maxminusdxminus;
  224. for(int i=0;i<N;i++){
  225. if (*ptr2maxminusdxminus < 0) {
  226. *ptr2maxminusdxminus = 0;
  227. } else {
  228. *ptr2maxminusdxminus *= *ptr2maxminusdxminus;
  229. }
  230. ptr2maxminusdxminus++;
  231. }
  233. float *mindxplus, *ptr2mindxplus;
  234. if((mindxplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_MDX+n");
  235. for(int i=0;i<N;i++)mindxplus[i]=dxplus[i];
  236. ptr2mindxplus = mindxplus;
  237. for(int i=0;i<N;i++){
  238. if (*ptr2mindxplus > 0) {
  239. *ptr2mindxplus = 0;
  240. } else {
  241. *ptr2mindxplus *= *ptr2mindxplus;
  242. }
  243. ptr2mindxplus++;
  244. }
  246. float *minminusdxminus, *ptr2minminusdxminus;
  247. if((minminusdxminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_-MDX-n");
  248. for(int i=0;i<N;i++)minminusdxminus[i]=-dxminus[i];
  249. ptr2minminusdxminus = minminusdxminus;
  250. for(int i=0;i<N;i++){
  251. if (*ptr2minminusdxminus > 0) {
  252. *ptr2minminusdxminus = 0;
  253. } else {
  254. *ptr2minminusdxminus *= *ptr2minminusdxminus;
  255. }
  256. ptr2minminusdxminus++;
  257. }
  259. //////// DY ////////
  261. float *dy, *ptr2dy;
  262. if((dy=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dyn");
  263. ptr2dy=dy;
  264. ptr2phi=phi;
  265. for(int c=0;c<imageW;c++){
  266. *ptr2dy = 0;
  267. ptr2dy++;
  268. ptr2phi++;
  269. }
  270. for(int r=0;r<imageH-2;r++){
  271. for(int c=0;c<imageW;c++){
  272. *ptr2dy = (*(ptr2phi-imageW) - *(ptr2phi+imageW))/2;
  273. ptr2dy++;
  274. ptr2phi++;
  275. }
  276. }
  277. for(int c=0;c<imageW;c++){
  278. *ptr2dy = 0;
  279. ptr2dy++;
  280. ptr2phi++;
  281. }
  283. float *dyplus, *ptr2dyplus;
  284. if((dyplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dy+n");
  285. for(int i=0;i<N;i++)dyplus[i]=phi[i];
  286. ptr2dyplus=dyplus+N-1;
  287. for(int r=0;r<imageH-1;r++){
  288. for(int c=0;c<imageW;c++){
  289. *ptr2dyplus = *(ptr2dyplus-imageW) - *ptr2dyplus;
  290. ptr2dyplus--;
  291. }
  292. }
  293. for(int c=0;c<imageW;c++){
  294. *ptr2dyplus = 0;
  295. ptr2dyplus--;
  296. }
  298. float *dyminus, *ptr2dyminus;
  299. if((dyminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dy+n");
  300. for(int i=0;i<N;i++)dyminus[i]=phi[i];
  301. ptr2dyminus=dyminus;
  302. for(int r=0;r<imageH-1;r++){
  303. for(int c=0;c<imageW;c++){
  304. *ptr2dyminus = *ptr2dyminus - *(ptr2dyminus+imageW);
  305. ptr2dyminus++;
  306. }
  307. }
  308. for(int c=0;c<imageW;c++){
  309. *ptr2dyminus = 0;
  310. ptr2dyminus++;
  311. }
  313. float *dyplusx, *ptr2dyplusx;
  314. if((dyplusx=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dyplusxn");
  315. ptr2dyplusx=dyplusx;
  316. ptr2phi=phi;
  317. for(int c=0;c<imageW;c++){
  318. *ptr2dyplusx = 0;
  319. ptr2dyplusx++;
  320. ptr2phi++;
  321. }
  322. for(int r=0;r<imageH-2;r++){
  323. for(int c=0;c<imageW-1;c++){
  324. ptr2phi++;
  325. *ptr2dyplusx = (*(ptr2phi-imageW) - *(ptr2phi+imageW))/2;
  326. ptr2dyplusx++;
  327. }
  328. *ptr2dyplusx=0;
  329. ptr2dyplusx++;
  330. ptr2phi++;
  331. }
  332. for(int c=0;c<imageW;c++){
  333. *ptr2dyplusx = 0;
  334. ptr2dyplusx++;
  335. }
  337. float *dyminusx, *ptr2dyminusx;
  338. if((dyminusx=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_dyminusxn");
  339. ptr2dyminusx=dyminusx;
  340. ptr2phi=phi;
  341. for(int c=0;c<imageW;c++){
  342. *ptr2dyminusx = 0;
  343. ptr2dyminusx++;
  344. ptr2phi++;
  345. }
  346. for(int r=0;r<imageH-2;r++){
  347. *ptr2dyminusx=0;
  348. for(int c=0;c<imageW-1;c++){
  349. ptr2dyminusx++;
  350. *ptr2dyminusx = (*(ptr2phi-imageW) - *(ptr2phi+imageW))/2;
  351. ptr2phi++;
  352. }
  353. ptr2dyminusx++;
  354. ptr2phi++;
  355. }
  356. for(int c=0;c<imageW;c++){
  357. *ptr2dyminusx = 0;
  358. ptr2dyminusx++;
  359. }
  361. float *maxdyplus, *ptr2maxdyplus;
  362. if((maxdyplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_Mdy+n");
  363. for(int i=0;i<N;i++)maxdyplus[i]=dyplus[i];
  364. ptr2maxdyplus = maxdyplus;
  365. for(int i=0;i<N;i++){
  366. if (*ptr2maxdyplus < 0) {
  367. *ptr2maxdyplus = 0;
  368. } else {
  369. *ptr2maxdyplus *= *ptr2maxdyplus;
  370. }
  371. ptr2maxdyplus++;
  372. }
  374. float *maxminusdyminus, *ptr2maxminusdyminus;
  375. if((maxminusdyminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_-Mdy-n");
  376. for(int i=0;i<N;i++)maxminusdyminus[i]=-dyminus[i];
  377. ptr2maxminusdyminus = maxminusdyminus;
  378. for(int i=0;i<N;i++){
  379. if (*ptr2maxminusdyminus < 0) {
  380. *ptr2maxminusdyminus = 0;
  381. } else {
  382. *ptr2maxminusdyminus *= *ptr2maxminusdyminus;
  383. }
  384. ptr2maxminusdyminus++;
  385. }
  387. float *mindyplus, *ptr2mindyplus;
  388. if((mindyplus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_Mdy+n");
  389. for(int i=0;i<N;i++)mindyplus[i]=dyplus[i];
  390. ptr2mindyplus = mindyplus;
  391. for(int i=0;i<N;i++){
  392. if (*ptr2mindyplus > 0) {
  393. *ptr2mindyplus = 0;
  394. } else {
  395. *ptr2mindyplus *= *ptr2mindyplus;
  396. }
  397. ptr2mindyplus++;
  398. }
  400. float *minminusdyminus, *ptr2minminusdyminus;
  401. if((minminusdyminus=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("ME_-Mdy-n");
  402. for(int i=0;i<N;i++)minminusdyminus[i]=-dyminus[i];
  403. ptr2minminusdyminus = minminusdyminus;
  404. for(int i=0;i<N;i++){
  405. if (*ptr2minminusdyminus > 0) {
  406. *ptr2minminusdyminus = 0;
  407. } else {
  408. *ptr2minminusdyminus *= *ptr2minminusdyminus;
  409. }
  410. ptr2minminusdyminus++;
  411. }
  413. float *gradphimax;
  414. if((gradphimax=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("GRADPHIMAXn");
  415. for(int i=0;i<N;i++){
  416. gradphimax[i]=sqrt((sqrt(maxdxplus[i]+maxminusdxminus[i]))*(sqrt(maxdxplus[i]+maxminusdxminus[i]))+(sqrt(maxdyplus[i]+maxminusdyminus[i]))*(sqrt(maxdyplus[i]+maxminusdyminus[i])));
  417. }
  419. float *gradphimin;
  420. if((gradphimin=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("GRADPHIminn");
  421. for(int i=0;i<N;i++){
  422. gradphimin[i]=sqrt((sqrt(mindxplus[i]+minminusdxminus[i]))*(sqrt(mindxplus[i]+minminusdxminus[i]))+(sqrt(mindyplus[i]+minminusdyminus[i]))*(sqrt(mindyplus[i]+minminusdyminus[i])));
  423. }
  425. float *nplusx;
  426. if((nplusx=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("nplusxn");
  427. for(int i=0;i<N;i++){
  428. nplusx[i]= dxplus[i] / sqrt(FLT_EPSILON + (dxplus[i]*dxplus[i]) + ((dyplusx[i] + dy[i])*(dyplusx[i] + dy[i])*0.25) );
  429. }
  430. float *nplusy;
  431. if((nplusy=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("nplusyn");
  432. for(int i=0;i<N;i++){
  433. nplusy[i]= dyplus[i] / sqrt(FLT_EPSILON + (dyplus[i]*dyplus[i]) + ((dxplusy[i] + dx[i])*(dxplusy[i] + dx[i])*0.25) );
  434. }
  436. float *nminusx;
  437. if((nminusx=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("nminusxn");
  438. for(int i=0;i<N;i++){
  439. nminusx[i]= dxminus[i] / sqrt(FLT_EPSILON + (dxminus[i]*dxminus[i]) + ((dyminusx[i] + dy[i])*(dyminusx[i] + dy[i])*0.25) );
  440. }
  442. float *nminusy;
  443. if((nminusy=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("nminusyn");
  444. for(int i=0;i<N;i++){
  445. nminusy[i]= dyminus[i] / sqrt(FLT_EPSILON + (dyminus[i]*dyminus[i]) + ((dxminusy[i] + dx[i])*(dxminusy[i] + dx[i])*0.25) );
  446. }
  448. float *curvature, *ptr2curvature;
  449. if((curvature=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("curvaturen");
  450. for(int i=0;i<N;i++){
  451. curvature[i]= ((nplusx[i]-nminusx[i])+(nplusy[i]-nminusy[i])/2);
  452. }
  455. float *F, *ptr2F, *ptr2D;
  456. if((F=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("Fn");
  457. ptr2F=F;
  458. ptr2D=D;
  459. ptr2curvature=curvature;
  460. for(int i=0;i<N;i++){
  461. *ptr2F = -(ALPHA * (*ptr2D)) + ( (1-ALPHA) * (*ptr2curvature));
  462. ptr2F++;
  463. ptr2D++;
  464. ptr2curvature++;
  465. }
  467. float *gradphi, *ptr2gradphi, *ptr2gradphimax, *ptr2gradphimin;
  468. if((gradphi=(float *)malloc(imageW*imageH*sizeof(float)))==NULL)printf("GRADPHIn");
  469. ptr2gradphi=gradphi;
  470. ptr2gradphimax=gradphimax;
  471. ptr2gradphimin=gradphimin;
  472. ptr2F=F;
  473. for(int i=0; i<N; i++){
  474. if(*ptr2F>0){
  475. *ptr2gradphi = *ptr2gradphimax;
  476. } else {
  477. *ptr2gradphi = *ptr2gradphimin;
  478. }
  479. ptr2gradphi++;
  480. ptr2F++;
  481. ptr2gradphimax++;
  482. ptr2gradphimin++;
  483. }
  485. ptr2phi=phi;
  486. ptr2gradphi=gradphi;
  487. ptr2F=F;
  488. for(int i=0; i<N; i++){
  489. *ptr2phi = *ptr2phi + DT * (*ptr2F) * (*ptr2gradphi);
  490. ptr2phi++;
  491. ptr2gradphi++;
  492. ptr2F++;
  493. }
  498. //printf("Freeing Memoryn");
  503. free(dx);
  504. free(dxplus);
  505. free(dxminus);
  506. free(dxminusy);
  507. free(maxdxplus);
  508. free(maxminusdxminus);
  509. free(mindxplus);
  510. free(minminusdxminus);
  511. //
  512. free(dy);
  513. free(dyplus);
  514. free(dyminus);
  515. free(dyminusx);
  516. free(maxdyplus);
  517. free(maxminusdyminus);
  518. free(mindyplus);
  519. free(minminusdyminus);
  520. //
  521. free(gradphi);
  522. free(gradphimax);
  523. free(gradphimin);
  524. //
  525. free(nplusx);
  526. free(nplusy);
  527. free(nminusx);
  528. free(nminusy);
  529. free(curvature);
  531. free(F);
  533. }
  538. int main(int argc, char** argv){
  540. const char *image_path = IMAGE;
  542. LoadBMPFile(&h_Src, &imageW, &imageH, image_path);
  543. D = (float *)malloc(imageW*imageH*sizeof(float));
  545. //printf("Input Imagen");
  546. for(int r=0;r<imageH;r++){
  547. for(int c=0;c<imageW;c++){
  548. D[r*imageW+c] = h_Src[r*imageW+c].x;
  549. /*printf("%3.0f ", D[r*imageW+c]);*/
  550. }
  551. //printf("n");
  552. }
  554. N = imageW*imageH;
  556. float *ptr2D;
  557. ptr2D=D;
  558. for(int i=0;i<N;i++){
  559. *ptr2D = EPSILON - abs(*ptr2D - THRESHOLD);
  560. ptr2D++;
  561. }
  563. init_phi();
  565. // Set up CUDA Timer
  566. cutCreateTimer(&Timer);
  567. cutStartTimer(Timer);
  569. for(its=0;its<ITERATIONS;its++){
  570. update_phi();
  571. if(its%50==0)printf("Iteration %3d Total Time: %3.2fn", its, 0.001*cutGetTimerValue(Timer));
  572. }
  573. }