API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jul 6th, 2015
206
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 15.39 KB | None | 0 0
raw download clone embed print report
  1. *fima = (float*)malloc(sizeof(float) * totalsize);
  2. means = (float*)malloc(sizeof(float) * totalsize);
  3. variances = (float*)malloc(sizeof(float) * totalsize);
  4. Estimate = (float*)malloc(sizeof(float) * totalsize);
  5. Label = (unsigned short*)malloc(sizeof(unsigned short) * totalsize);
  6.  
  7. /* Pierrick Coupe - pierrick.coupe@gmail.com */
  8. /* Jose V. Manjon - jmanjon@fis.upv.es */
  9. /* Brain Imaging Center, Montreal Neurological Institute. */
  10. /* Mc Gill University */
  11. /* */
  12. /* Copyright (C) 2010 Pierrick Coupe and Jose V. Manjon */
  13.  
  14. /* Details on ONLM filter */
  15. /***************************************************************************
  16. * The ONLM filter is described in: *
  17. * *
  18. * P. Coupe, P. Yger, S. Prima, P. Hellier, C. Kervrann, C. Barillot. *
  19. * An Optimized Blockwise Non Local Means Denoising Filter for 3D Magnetic*
  20. * Resonance Images. IEEE Transactions on Medical Imaging, 27(4):425-441, *
  21. * Avril 2008 *
  22. ***************************************************************************/
  23. /***************************************************************************
  24. * This file was modified by Haider Khan - haiderriazkhan@hotmail.com *
  25. * *
  26. * The ONLM C program can now be called from Java via JNA *
  27. ***************************************************************************/
  28.  
  29.  
  30. #include "math.h"
  31. #include <stdlib.h>
  32. #include <string.h>
  33. #include <stdio.h>
  34.  
  35.  
  36.  
  37. /* Multithreading stuff*/
  38. #ifdef _WIN32
  39. #include <windows.h>
  40. #include <process.h>
  41. #else
  42. #include <pthread.h>
  43. #endif
  44.  
  45. #include <stdbool.h>
  46.  
  47.  
  48.  
  49.  
  50. typedef struct{
  51. int rows;
  52. int cols;
  53. int slices;
  54. float * in_image;
  55. float * ref_image;
  56. float * means_image;
  57. float * var_image;
  58. float * estimate;
  59. unsigned short * label;
  60. int ini;
  61. int fin;
  62. int radioB;
  63. int radioS;
  64. float *sigma;
  65. float beta;
  66. }myargument;
  67.  
  68.  
  69. float max;
  70.  
  71. /* Function which compute the weighted average for one block */
  72. void Average_block(float *ima,int x,int y,int z,int neighborhoodsize,float *average, float weight, int sx,int sy,int sz)
  73. {
  74. int x_pos,y_pos,z_pos;
  75. bool is_outside;
  76. int a,b,c,ns,sxy,count;
  77.  
  78. ns=2*neighborhoodsize+1;
  79. sxy=sx*sy;
  80.  
  81. count = 0;
  82.  
  83. for (c = 0; c<ns;c++)
  84. {
  85. z_pos = z+c-neighborhoodsize;
  86. for (b = 0; b<ns;b++)
  87. {
  88. y_pos = y+b-neighborhoodsize;
  89. for (a = 0; a<ns;a++)
  90. {
  91. x_pos = x+a-neighborhoodsize;
  92.  
  93. is_outside = false;
  94. if ((x_pos < 0) || (x_pos > sx-1)) is_outside = true;
  95. if ((y_pos < 0) || (y_pos > sy-1)) is_outside = true;
  96. if ((z_pos < 0) || (z_pos > sz-1)) is_outside = true;
  97.  
  98.  
  99. if (is_outside)
  100. average[count] = average[count] + ima[z*(sxy)+(y*sx)+x]*weight;
  101. else
  102. average[count] = average[count] + ima[z_pos*(sxy)+(y_pos*sx)+x_pos]*weight;
  103.  
  104. count++;
  105. }
  106. }
  107. }
  108. }
  109.  
  110. /* Function which computes the value assigned to each voxel */
  111. void Value_block(float *Estimate, unsigned short *Label,int x,int y,int z,int neighborhoodsize,float *average, float global_sum, int sx,int sy,int sz)
  112. {
  113. int x_pos,y_pos,z_pos;
  114. int ret;
  115. bool is_outside;
  116. float value = 0.0;
  117. unsigned short label = 0;
  118. int count=0 ;
  119. int a,b,c,ns,sxy;
  120.  
  121. extern bool rician;
  122.  
  123. ns=2*neighborhoodsize+1;
  124. sxy=sx*sy;
  125.  
  126. for (c = 0; c<ns;c++)
  127. {
  128. z_pos = z+c-neighborhoodsize;
  129. for (b = 0; b<ns;b++)
  130. {
  131. y_pos = y+b-neighborhoodsize;
  132. for (a = 0; a<ns;a++)
  133. {
  134. x_pos = x+a-neighborhoodsize;
  135.  
  136. is_outside = false;
  137. if ((x_pos < 0) || (x_pos > sx-1)) is_outside = true;
  138. if ((y_pos < 0) || (y_pos > sy-1)) is_outside = true;
  139. if ((z_pos < 0) || (z_pos > sz-1)) is_outside = true;
  140.  
  141. if (!is_outside)
  142. {
  143. value = Estimate[z_pos*(sxy)+(y_pos*sx)+x_pos];
  144.  
  145. value = value + (average[count]/global_sum);
  146.  
  147. label = Label[(x_pos + y_pos*sx + z_pos*sxy)];
  148. Estimate[z_pos*(sxy)+(y_pos*sx)+x_pos] = value;
  149. Label[(x_pos + y_pos*sx + z_pos *sxy)] = label +1;
  150. }
  151. count++;
  152. }
  153. }
  154. }
  155. }
  156.  
  157. float distance(float* ima,int x,int y,int z,int nx,int ny,int nz,int f,int sx,int sy,int sz)
  158. {
  159. float d,acu,distancetotal;
  160. int i,j,k,ni1,nj1,ni2,nj2,nk1,nk2;
  161.  
  162. distancetotal=0;
  163.  
  164. for(k=-f;k<=f;k++)
  165. {
  166. nk1=z+k;
  167. nk2=nz+k;
  168. if(nk1<0) nk1=-nk1;
  169. if(nk2<0) nk2=-nk2;
  170. if(nk1>=sz) nk1=2*sz-nk1-1;
  171. if(nk2>=sz) nk2=2*sz-nk2-1;
  172.  
  173. for(j=-f;j<=f;j++)
  174. {
  175. nj1=y+j;
  176. nj2=ny+j;
  177. if(nj1<0) nj1=-nj1;
  178. if(nj2<0) nj2=-nj2;
  179. if(nj1>=sy) nj1=2*sy-nj1-1;
  180. if(nj2>=sy) nj2=2*sy-nj2-1;
  181.  
  182. for(i=-f;i<=f;i++)
  183. {
  184. ni1=x+i;
  185. ni2=nx+i;
  186. if(ni1<0) ni1=-ni1;
  187. if(ni2<0) ni2=-ni2;
  188. if(ni1>=sx) ni1=2*sx-ni1-1;
  189. if(ni2>=sx) ni2=2*sx-ni2-1;
  190.  
  191. distancetotal = distancetotal + ((ima[nk1*(sx*sy)+(nj1*sx)+ni1]-ima[nk2*(sx*sy)+(nj2*sx)+ni2])*(ima[nk1*(sx*sy)+(nj1*sx)+ni1]-ima[nk2*(sx*sy)+(nj2*sx)+ni2]));
  192. }
  193. }
  194. }
  195.  
  196. acu=(2*f+1)*(2*f+1)*(2*f+1);
  197. d=distancetotal/acu;
  198.  
  199. return d;
  200.  
  201. }
  202.  
  203. void* ThreadFunc( void* pArguments )
  204. {
  205. float *Estimate,*ima,*means,*variances,*sigma,*ref,*average,epsilon,totalweight,wmax,t1,d,w;
  206. float beta;
  207. unsigned short *Label;
  208. int rows,cols,slices,ini,fin,v,f,i,j,k,rc,ii,jj,kk,ni,nj,nk,Ndims;
  209.  
  210. myargument arg;
  211. arg=*(myargument *) pArguments;
  212.  
  213. rows=arg.rows;
  214. cols=arg.cols;
  215. slices=arg.slices;
  216. ini=arg.ini;
  217. fin=arg.fin;
  218. ima=arg.in_image;
  219. ref=arg.ref_image;
  220. means=arg.means_image;
  221. variances=arg.var_image;
  222. Estimate=arg.estimate;
  223. Label=arg.label;
  224. v=arg.radioB;
  225. f=arg.radioS;
  226. sigma=arg.sigma;
  227. beta= arg.beta;
  228.  
  229. epsilon = 0.0000001;
  230.  
  231. rc=rows*cols;
  232.  
  233. Ndims = (2*f+1)*(2*f+1)*(2*f+1);
  234.  
  235. average=(float*)malloc(Ndims*sizeof(float));
  236.  
  237. for(k=ini;k<fin;k+=2)
  238. {
  239. for(j=0;j<rows;j+=2)
  240. {
  241. for(i=0;i<cols;i+=2)
  242. {
  243. memset(average,0.0, sizeof(float) * Ndims );
  244.  
  245. totalweight=0.0;
  246.  
  247. wmax=0.0;
  248.  
  249.  
  250. for(kk=-v;kk<=v;kk++)
  251. {
  252. for(jj=-v;jj<=v;jj++)
  253. {
  254. for(ii=-v;ii<=v;ii++)
  255. {
  256. ni=i+ii;
  257. nj=j+jj;
  258. nk=k+kk;
  259.  
  260. if(ii==0 && jj==0 && kk==0) continue;
  261.  
  262. if(ni>=0 && nj>=0 && nk>=0 && ni<cols && nj<rows && nk<slices)
  263. {
  264. t1 = ((2 * means[k*rc+(j*cols)+i] * means[nk*rc+(nj*cols)+ni]+ epsilon) / ( means[k*rc+(j*cols)+i]*means[k*rc+(j*cols)+i] + means[nk*rc+(nj*cols)+ni]*means[nk*rc+(nj*cols)+ni] + epsilon) ) * ((2 * sqrt(variances[k*rc+(j*cols)+i]) * sqrt(variances[nk*rc+(nj*cols)+ni]) + epsilon) / (variances[k*rc+(j*cols)+i] + variances[nk*rc+(nj*cols)+ni] + epsilon));
  265. if(t1 > 0.9)
  266. {
  267.  
  268. d=distance(ref,i,j,k,ni,nj,nk,f,cols,rows,slices);
  269. if (d<=0) d = epsilon;
  270.  
  271. w = exp(-d/(beta*sigma[k*(rc)+(j*cols)+i]*sigma[k*(rc)+(j*cols)+i]));
  272.  
  273. if(w>wmax) wmax = w;
  274.  
  275.  
  276. Average_block(ima,ni,nj,nk,f,average,w,cols,rows,slices);
  277. totalweight = totalweight + w;
  278.  
  279. }
  280.  
  281.  
  282. }
  283. }
  284. }
  285.  
  286. }
  287.  
  288. if(wmax==0.0) wmax=1.0;
  289.  
  290. Average_block(ima,i,j,k,f,average,wmax,cols,rows,slices);
  291.  
  292. totalweight = totalweight + wmax;
  293.  
  294.  
  295. if(totalweight != 0.0)
  296. Value_block(Estimate,Label,i,j,k,f,average,totalweight,cols,rows,slices);
  297.  
  298. }
  299. }
  300. }
  301.  
  302. free(average);
  303.  
  304.  
  305. #ifdef _WIN32
  306. _endthreadex(0);
  307. return 0;
  308. #else
  309. pthread_exit(0);
  310. return 0;
  311. #endif
  312.  
  313. }
  314.  
  315.  
  316.  
  317. void cleanup(float* pVals)
  318. {
  319. free(pVals);
  320. }
  321.  
  322.  
  323.  
  324. void ONLM(float* ima , int v , int f , float* sigma , float beta , float* ref , int rows, int cols, int slices, float** fima)
  325.  
  326. {
  327.  
  328.  
  329. float *means, *variances, *Estimate;
  330. unsigned short *Label;
  331. float w,totalweight,wmax,d,mean,var,t1,t2,epsilon,label,estimate;
  332.  
  333. int i,j,k,ii,jj,kk,ni,nj,nk,indice,Nthreads,ini,fin;
  334.  
  335.  
  336.  
  337. myargument *ThreadArgs;
  338.  
  339. #ifdef _WIN32
  340. HANDLE *ThreadList; // Handles to the worker threads
  341. #else
  342. pthread_t * ThreadList;
  343. #endif
  344.  
  345.  
  346.  
  347.  
  348. const int dims[3] = {cols, rows, slices};
  349.  
  350.  
  351. int totalsize = (rows*cols*slices);
  352. *fima = (float*)malloc(sizeof(float) * totalsize);
  353. means = (float*)malloc(sizeof(float) * totalsize);
  354. variances = (float*)malloc(sizeof(float) * totalsize);
  355. Estimate = (float*)malloc(sizeof(float) * totalsize);
  356. Label = (unsigned short*)malloc(sizeof(unsigned short) * totalsize);
  357.  
  358.  
  359.  
  360.  
  361. for (i = 0; i < dims[2] *dims[1] * dims[0];i++)
  362. {
  363. Estimate[i] = 0.0;
  364. Label[i] = 0.0;
  365. (*fima)[i] = 0.0;
  366. }
  367.  
  368.  
  369. max=0;
  370. for(k=0;k<dims[2];k++)
  371. {
  372. for(i=0;i<dims[1];i++)
  373. {
  374. for(j=0;j<dims[0];j++)
  375. {
  376. mean=0;
  377. indice=0;
  378. for(ii=-1;ii<=1;ii++)
  379. {
  380. for(jj=-1;jj<=1;jj++)
  381. {
  382. for(kk=-1;kk<=1;kk++)
  383. {
  384. ni=i+ii;
  385. nj=j+jj;
  386. nk=k+kk;
  387.  
  388. if(ni<0) ni=-ni;
  389. if(nj<0) nj=-nj;
  390. if(nk<0) nk=-nk;
  391. if(ni>=dims[1]) ni=2*dims[1]-ni-1;
  392. if(nj>=dims[0]) nj=2*dims[0]-nj-1;
  393. if(nk>=dims[2]) nk=2*dims[2]-nk-1;
  394.  
  395.  
  396. mean = mean + ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj];
  397. indice=indice+1;
  398.  
  399. }
  400. }
  401. }
  402. mean=mean/indice;
  403. means[k*(dims[0]*dims[1])+(i*dims[0])+j]=mean;
  404. if (mean > max) max =mean;
  405. }
  406. }
  407. }
  408.  
  409.  
  410. for(k=0;k<dims[2];k++)
  411. {
  412. for(i=0;i<dims[1];i++)
  413. {
  414. for(j=0;j<dims[0];j++)
  415. {
  416. var=0;
  417. indice=0;
  418. for(ii=-1;ii<=1;ii++)
  419. {
  420. for(jj=-1;jj<=1;jj++)
  421. {
  422. for(kk=-1;kk<=1;kk++)
  423. {
  424. ni=i+ii;
  425. nj=j+jj;
  426. nk=k+kk;
  427. if(ni>=0 && nj>=0 && nk>0 && ni<dims[1] && nj<dims[0] && nk<dims[2])
  428. {
  429. var = var + (ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj]-means[k*(dims[0]*dims[1])+(i*dims[0])+j])*(ref[nk*(dims[0]*dims[1])+(ni*dims[0])+nj]-means[k*(dims[0]*dims[1])+(i*dims[0])+j]);
  430. indice=indice+1;
  431. }
  432. }
  433. }
  434. }
  435. var=var/(indice-1);
  436. variances[k*(dims[0]*dims[1])+(i*dims[0])+j]=var;
  437. }
  438. }
  439. }
  440.  
  441. Nthreads=8;
  442.  
  443.  
  444. #ifdef _WIN32
  445.  
  446. // Reserve room for handles of threads in ThreadList
  447. ThreadList = (HANDLE*)malloc(Nthreads* sizeof( HANDLE ));
  448. ThreadArgs = (myargument*) malloc( Nthreads*sizeof(myargument));
  449.  
  450. for (i=0; i<Nthreads; i++)
  451. {
  452. /* Make Thread Structure*/
  453. ini=(i*dims[2])/Nthreads;
  454. fin=((i+1)*dims[2])/Nthreads;
  455. ThreadArgs[i].cols=dims[0];
  456. ThreadArgs[i].rows=dims[1];
  457. ThreadArgs[i].slices=dims[2];
  458. ThreadArgs[i].in_image=ima;
  459. ThreadArgs[i].ref_image=ref;
  460. ThreadArgs[i].var_image=variances;
  461. ThreadArgs[i].means_image=means;
  462. ThreadArgs[i].estimate=Estimate;
  463. ThreadArgs[i].label=Label;
  464. ThreadArgs[i].ini=ini;
  465. ThreadArgs[i].fin=fin;
  466. ThreadArgs[i].radioB=v;
  467. ThreadArgs[i].radioS=f;
  468. ThreadArgs[i].sigma=sigma;
  469. ThreadArgs[i].beta=beta;
  470. ThreadList[i] = (HANDLE)_beginthreadex( NULL, 0, &ThreadFunc, &ThreadArgs[i] , 0, NULL );
  471. }
  472.  
  473. for (i=0; i<Nthreads; i++) { WaitForSingleObject(ThreadList[i], INFINITE); }
  474. for (i=0; i<Nthreads; i++) { CloseHandle( ThreadList[i] ); }
  475.  
  476. #else
  477.  
  478. /* Reserve room for handles of threads in ThreadList*/
  479. ThreadList = (pthread_t *) calloc(Nthreads,sizeof(pthread_t));
  480. ThreadArgs = (myargument*) calloc( Nthreads,sizeof(myargument));
  481.  
  482. for (i=0; i<Nthreads; i++)
  483. {
  484. /* Make Thread Structure*/
  485. ini=(i*dims[2])/Nthreads;
  486. fin=((i+1)*dims[2])/Nthreads;
  487. ThreadArgs[i].cols=dims[0];
  488. ThreadArgs[i].rows=dims[1];
  489. ThreadArgs[i].slices=dims[2];
  490. ThreadArgs[i].in_image=ima;
  491. ThreadArgs[i].ref_image=ref;
  492. ThreadArgs[i].var_image=variances;
  493. ThreadArgs[i].means_image=means;
  494. ThreadArgs[i].estimate=Estimate;
  495. ThreadArgs[i].label=Label;
  496. ThreadArgs[i].ini=ini;
  497. ThreadArgs[i].fin=fin;
  498. ThreadArgs[i].radioB=v;
  499. ThreadArgs[i].radioS=f;
  500. ThreadArgs[i].sigma=sigma;
  501. ThreadArgs[i].beta=beta;
  502. }
  503.  
  504. for (i=0; i<Nthreads; i++)
  505. {
  506.  
  507. if(pthread_create(&ThreadList[i], NULL, ThreadFunc,&ThreadArgs[i]))
  508. {
  509. printf("Threads cannot be createdn");
  510. exit(1);
  511. }
  512.  
  513. }
  514.  
  515.  
  516.  
  517. for (i=0; i<Nthreads; i++)
  518. {
  519. pthread_join(ThreadList[i],NULL);
  520. }
  521.  
  522. #endif
  523.  
  524. free(ThreadArgs);
  525. free(ThreadList);
  526.  
  527. free(variances);
  528. free(means);
  529. free(sigma);
  530.  
  531.  
  532. label = 0.0;
  533. estimate = 0.0;
  534.  
  535.  
  536. /* Aggregation of the estimators (i.e. means computation) */
  537. for (k = 0; k < dims[2]; k++ )
  538. {
  539. for (i = 0; i < dims[1]; i++ )
  540. {
  541. for (j = 0; j < dims[0]; j++ )
  542. {
  543. label = Label[k*(dims[0]*dims[1])+(i*dims[0])+j];
  544. if (label == 0.0)
  545. {
  546. // Corrected
  547. (*fima)[k*(dims[0]*dims[1])+(i*dims[0])+j] = ima[k*(dims[0]*dims[1])+(i*dims[0])+j];
  548.  
  549. }
  550. else
  551. {
  552. estimate = Estimate[k*(dims[0]*dims[1])+(i*dims[0])+j];
  553. estimate = (estimate/label);
  554. (*fima)[k*(dims[0]*dims[1])+(i*dims[0])+j]=estimate;
  555.  
  556. }
  557. }
  558. }
  559. }
  560.  
  561.  
  562.  
  563.  
  564.  
  565. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!