Fixed entghost bug

entghost.c

@@ -190,14 +190,16 @@ void phase_correct(float *in, float *out, int ppl, int lpi,int nseg,float ph1)
 double calc_entropy(float *data,long imsize)
 {
   long l;
-  double s;
+  double s,t;
   int hist[100];
 
   calc_hist(data,imsize,egr_min,egr_max,100,hist);
   
   s=0;
   for(l=0;l<100;l++){
-    s-=((double)hist[l]/imsize)*log((double)hist[l]/imsize);
+    t=((double)hist[l]/(double)imsize)*log((double)hist[l]/(double)imsize);
+    if(!isnan(t))
+      s-=t;
   }
   return(s/log(100.0)); 
 }
@@ -330,12 +332,12 @@ void iepi_entropy_ghost_reduction(float *data,int ppl,int lpi,int ipv,int nseg)
   egr_ppl = ppl;
   egr_lpi = lpi;
   egr_seg = nseg;
-  egr_tmp = (float *)malloc(ppl*lpi*2*sizeof(float));
 
   ColumnVector result(4*nseg);
   ColumnVector tol(4*nseg);
   ColumnVector bound(4*nseg);
 
+  result=0;
   tol=0.01;
   bound = 10;
 
@@ -344,14 +346,15 @@ void iepi_entropy_ghost_reduction(float *data,int ppl,int lpi,int ipv,int nseg)
     egr_data = &data[i*ppl*lpi*2];
     result=0;
     egr_calc = 0;
-    optimise(result,nseg*3,tol,iepi_cost_function,iter,100,bound);
+    optimise(result,nseg*4,tol,iepi_cost_function,iter,100,bound);
     iepi_rephase(egr_data,egr_data,ppl,lpi,nseg,result);
+    //cout << result.t();
   }
-  free(egr_tmp);
 }
 
 float iepi_cost_function(const ColumnVector& input)
 {
+  egr_tmp = (float *)calloc(egr_ppl*egr_lpi*2,sizeof(float));
   iepi_rephase(egr_data, egr_tmp, egr_ppl, egr_lpi, egr_seg, input);
   image_fftw_c(egr_tmp, egr_ppl, egr_lpi, 1, 1);
   lv_modulus(egr_tmp,egr_ppl*egr_lpi*2);
@@ -362,6 +365,8 @@ float iepi_cost_function(const ColumnVector& input)
     egr_calc = 1;
   }
   float e = calc_entropy(egr_tmp,egr_ppl*egr_lpi);
+  free(egr_tmp);
+  
   return(e);
 }
 

process.c

@@ -592,7 +592,7 @@ int process(char *infile,char *outfile,struct flags *proc)
   /*     S T A R T   O F    M A I N   P R O C E S S I N G   L O O P       */
   /* ==================================================================== */
   /* ==================================================================== */
- 
+
   /* Loop over volumes */
   maxset=0;
   for(i=proc->start;i<proc->start+proc->num;i++){