src/Tools/ctfInfo/ctfDeterminationByBayes/src/ctfDeterminationByBayes.c

   1 /*
   2 # ctfDeterminationByBayes : $Revision$
   3 # $Date$
   4 # Created by $Author$
   5 # Usage : ctfDeterminationByBayes
   6 # Attention
   7 #   $Loccker$
   8 #       $State$
   9 #
  10 */
  11 #include <stdio.h>
  12 #include <stdlib.h>
  13 #include <string.h>
  14 #include <math.h>
  15 #define GLOBAL_DECLARATION
  16 #include "../inc/config.h"
  17 #include "eosBayes.h"
  18 #include "ctfInfo.h"
  19 #include "ctffuncforbayes.h"
  20
  21
  22 #define DEBUG
  23 #include "genUtil.h"
  24 #include "ctfDeterminationByBayes.h"
  25 #include "String.h"
  26
  27 /*
  28 Example:
  29 typedef struct lctfDeterminationByBayesInfo {
  30         float a;
  31         int   b;
  32 } lctfDeterminationByBayesInfo;
  33
  34 typedef enum lctfDeterminationByBayesMode {
  35         a=0,
  36         b=1
  37 } lctfDeterminationByBayesMode;
  38 */
  39
  40 void DefocusEstimationLoop(ctfDeterminationByBayesInfo info, eosBayes* out, ctfInfo* ctfinfo, ctfforbayes* ctf, intensityforbayes* intensity, defocuslikelihoodforbayes* likelihood, int count, int mode);
  41
  42 int
  43 main(int argc, char* argv[])
  44 {
  45
  46   int i;
  47
  48   ctfDeterminationByBayesInfo info;
  49   ctfInfo ctfinfo;
  50   eosBayesDistributionInfo llinfo;
  51
  52   eosBayes out;
  53
  54   ctfforbayes ctf;
  55   intensityforbayes intensity;
  56   defocuslikelihoodforbayes likelihood;
  57
  58   float Rmax;
  59   float Rmin;
  60   float Rintermediate;
  61   float Rfirstpeak;
  62
  63
  64   init0(&info);
  65   argCheck(&info, argc, argv);
  66   init1(&info);
  67
  68   DEBUGPRINT("Program Start\n");
  69   DEBUGPRINT("local binary:MTF plus\n");
  70
  71   ctfinfo.mode = 33;
  72
  73   /*For Prior*/
  74   leosBayesInit(&out, info.numDistribution, info.numLikelihood, info.mode);
  75   eosBayesCondition(info.fptIn, &out, &llinfo, 0);
  76   eosBayesPriorWrite(&out,0);
  77   DistributionChangeCheck(&out, 0, 0);
  78
  79
  80   /*For CTF and intensity*/
  81   ctfforbayesInit(info.fptIn2,&ctf,0);
  82   intensityforbayesInit(&intensity,ctf.n,info.MTFmode,info.Ainmode,0);
  83   spatialfrequencyRead(info.fptIn2,&ctf,info.rmax,info.rmin,0);
  84   ctfInfoRead(info.fptIn3,&ctfinfo,"",0);
  85   DEBUGPRINT6("kV:%f  Cs:%f  Ain:%f Cc:%f MTF:%f ctfMode:%ld\n",ctfinfo.kV,ctfinfo.Cs,ctfinfo.Ain,ctfinfo.Cc,ctfinfo.BofMTF,ctfinfo.mode);
  86
  87
  88   /*ctfInfoflagcheck*/
  89   DEBUGPRINT2("flagAliasing:%ld flagSampling:%ld\n",ctfinfo.flagAliasing,ctfinfo.flagSampling);
  90   DEBUGPRINT1("NyquistFrequency:%f\n",ctfinfo.NyquistFrequency);
  91
  92   ctfinfo.flagMolecEnvTable = 0;
  93   DEBUGPRINT1("flagMolecEnvTable:%ld\n",ctfinfo.flagMolecEnvTable);
  94   ctfinfo.flagElastic = 0;
  95   DEBUGPRINT1("flagElastic:%ld\n",ctfinfo.flagElastic);
  96
  97   ctfinfo.flagWithInElasticTable = 0;
  98   DEBUGPRINT1("flagWithInElasticEnvTable:%ld\n",ctfinfo.flagWithInElasticTable);
  99   ctfinfo.flagInElastic = 0;
 100   DEBUGPRINT1("flagInElastic:%ld\n",ctfinfo.flagInElastic);
 101   ctfinfo.flagWithInElastic = 0;
 102   DEBUGPRINT1("flagWithInElastic:%ld\n",ctfinfo.flagWithInElastic);
 103
 104   ctfinfo.flagVibration = 0;
 105   DEBUGPRINT1("flagVibration:%ld\n",ctfinfo.flagVibration);
 106   ctfinfo.VibrationMode = 0;
 107   DEBUGPRINT1("VibrationMode:%ld\n",ctfinfo.VibrationMode);
 108
 109   ctfinfo.Magnification = 0;
 110   DEBUGPRINT1("Magnification:%f\n",ctfinfo.Magnification);
 111   ctfinfo.flagMagnification = 0;
 112   DEBUGPRINT1("flagMagnification:%ld\n",ctfinfo.flagMagnification);
 113
 114
 115   DEBUGPRINT1("flagIn4:%ld\n",info.flagIn4);
 116
 117   Rmax = ctf.rmax;
 118   Rmin = ctf.rmin;
 119   Rintermediate = info.rintermediate;
 120   Rfirstpeak = info.rfirstpeak;
 121
 122   /*global carse search B*/
 123   DEBUGPRINT("carse search start\n");
 124   ctf.rmax = Rmax;
 125   ctf.rmin = Rintermediate;
 126   DEBUGPRINT2("set R range:%f~%f\n",ctf.rmin,ctf.rmax);
 127   DefocusEstimationLoop(info,&out,&ctfinfo,&ctf,&intensity,&likelihood,0,0);
 128   DEBUGPRINT("carse search finish\n");
 129   ctfforbayescheckIntensityOfMaxPosterior(ctf,out.posterior,out.numDistribution,&ctfinfo,intensity.MTFmode,intensity.Ainmode,1,1);
 130
 131   /*MTF serch start\n*/
 132   DEBUGPRINT("MTF,B search start\n");
 133   ctf.rmax = Rmax;
 134   ctf.rmin = Rintermediate;
 135   DEBUGPRINT2("set R range:%f~%f\n",ctf.rmin,ctf.rmax);
 136   DEBUGPRINT("----maxflagcheck----\n");
 137   #ifdef DEBUG
 138   for(i=0;i<out.numDistribution;i++){
 139         printf("dist%d:%d\n",i,(out.posterior[i].maxposteriorflag));
 140   }
 141   #endif
 142
 143   DEBUGPRINT("set canstantflag of Ain");
 144   out.posterior[0].constantflag = 1;
 145   out.posterior[1].constantflag = 1;
 146   out.posterior[3].constantflag = 1;
 147   out.posterior[5].constantflag = 1;
 148
 149   DEBUGPRINT("set resetflag of \n");
 150   out.prior[4].resetflag = 1;
 151   eosBayesProbabilityReset(info.fptIn,&out,10,0);
 152
 153   DEBUGPRINT("set rangechangeflag of B\n");
 154   eosBayesRangeChangeflagSet(&out,2,0.95,100);
 155   eosBayesProbabilityRangeChange(&out,0);
 156
 157   DistributionChangeCheck(&out, 1, 0);
 158
 159   DefocusEstimationLoop(info,&out,&ctfinfo,&ctf,&intensity,&likelihood,1,1);
 160   DEBUGPRINT("MTF,B search finish\n");
 161
 162   DEBUGPRINT("----maxflagcheck----\n");
 163   #ifdef DEBUG
 164   for(i=0;i<out.numDistribution;i++){
 165         printf("dist%d:%d\n",i,(out.posterior[i].maxposteriorflag));
 166   }
 167   #endif
 168
 169   ctfforbayescheckIntensityOfMaxPosterior(ctf,out.posterior,out.numDistribution,&ctfinfo,intensity.MTFmode,intensity.Ainmode,1,1);
 170
 171
 172   /*df,A B serch*/
 173   DEBUGPRINT("df, A, B search start\n");
 174   ctf.rmax = Rintermediate;
 175   ctf.rmin = Rmin;
 176   DEBUGPRINT2("set R range:%f~%f\n",ctf.rmin,ctf.rmax);
 177
 178   DEBUGPRINT("set maxposteriorflag of MTF,B\n");
 179   eosBayesMaxposteriorflagSet(&out,4);
 180
 181   DEBUGPRINT("set canstantflag of Ain");
 182   out.posterior[3].constantflag = 1;
 183   out.posterior[5].constantflag = 1;
 184
 185   DEBUGPRINT("set resetflag of df,k\n");
 186   out.prior[0].resetflag = 1;
 187   out.prior[1].resetflag = 1;
 188   eosBayesProbabilityReset(info.fptIn,&out,1,0);
 189
 190   DEBUGPRINT("set rangechangeflag of B\n");
 191   eosBayesRangeChangeflagSet(&out,2,0.95,20);
 192   eosBayesProbabilityRangeChange(&out,0);
 193
 194
 195   DistributionChangeCheck(&out, 4, 0);
 196   DEBUGPRINT("----maxflagcheck----\n");
 197   #ifdef DEBUG
 198   for(i=0;i<out.numDistribution;i++){
 199         printf("dist%d:%d\n",i,(out.posterior[i].maxposteriorflag));
 200   }
 201   #endif
 202
 203   DefocusEstimationLoop(info,&out,&ctfinfo,&ctf,&intensity,&likelihood,4,1);
 204   DEBUGPRINT("A,B,df search finish\n");
 205   ctfforbayescheckIntensityOfMaxPosterior(ctf,out.posterior,out.numDistribution,&ctfinfo,intensity.MTFmode,intensity.Ainmode,1,1);
 206
 207
 208   /*df,A, B 2 serch*/
 209   DEBUGPRINT("df, A, B search 2 start\n");
 210   ctf.rmax = Rintermediate;
 211   ctf.rmin = Rmin;
 212   DEBUGPRINT2("set R range:%f~%f\n",ctf.rmin,ctf.rmax);
 213
 214   DEBUGPRINT("set maxposteriorflag of MTF,B\n");
 215   eosBayesMaxposteriorflagSet(&out,4);
 216
 217   DEBUGPRINT("set canstantflag of Ain");
 218   out.posterior[3].constantflag = 1;
 219   out.posterior[5].constantflag = 1;
 220
 221   DEBUGPRINT("set resetflag of df,k\n");
 222   out.prior[0].resetflag = 1;
 223   eosBayesProbabilityReset(info.fptIn,&out,5,0);
 224
 225
 226   DEBUGPRINT("set rangechangeflag of A\n");
 227   eosBayesRangeChangeflagSet(&out,1,0.95,20);
 228   eosBayesRangeChangeflagSet(&out,2,0.95,20);
 229   eosBayesProbabilityRangeChange(&out,0);
 230
 231   DistributionChangeCheck(&out, 5, 0);
 232   DEBUGPRINT("----maxflagcheck----\n");
 233   #ifdef DEBUG
 234   for(i=0;i<out.numDistribution;i++){
 235         printf("dist%d:%d\n",i,(out.posterior[i].maxposteriorflag));
 236         printf("rangechange%d:%d\n",i,(out.posterior[i].rangechangeflag));
 237
 238   }
 239   #endif
 240
 241   DefocusEstimationLoop(info,&out,&ctfinfo,&ctf,&intensity,&likelihood,5,1);
 242   DEBUGPRINT("A,df,k search finish\n");
 243   ctfforbayescheckIntensityOfMaxPosterior(ctf,out.posterior,out.numDistribution,&ctfinfo,intensity.MTFmode,intensity.Ainmode,1,1);
 244
 245
 246   /*Defocus ref*/
 247   DEBUGPRINT("defocus refinement start\n");
 248   ctf.rmax = Rmax;
 249   ctf.rmin = Rfirstpeak;
 250   DEBUGPRINT2("set R range:%f~%f\n",ctf.rmin,ctf.rmax);
 251
 252   DEBUGPRINT("set maxposteriorflag of defocus,A,B,MTF\n");
 253   out.posterior[0].maxposteriorflag = 0;
 254   out.posterior[1].maxposteriorflag = 0;
 255   out.posterior[3].maxposteriorflag = 0;
 256   out.posterior[4].maxposteriorflag = 0;
 257   out.posterior[5].maxposteriorflag = 0;
 258   // eosBayesMaxposteriorflagSet(&out,4);
 259
 260   DEBUGPRINT("set canstantflag of Ain");
 261   out.posterior[3].constantflag = 0;
 262   out.posterior[5].constantflag = 0;
 263
 264   DEBUGPRINT("set resetflag of df,k\n");
 265   out.prior[3].resetflag = 1;
 266   out.prior[5].resetflag = 1;
 267   eosBayesProbabilityReset(info.fptIn,&out,1,0);
 268
 269   DEBUGPRINT("set rangechangeflag of A\n");
 270   eosBayesRangeChangeflagSet(&out,0,0.70,10);
 271   eosBayesRangeChangeflagSet(&out,1,0.95,20);
 272   eosBayesRangeChangeflagSet(&out,2,0.95,20);
 273   eosBayesRangeChangeflagSet(&out,4,0.40,5);
 274   eosBayesProbabilityRangeChange(&out,0);
 275
 276   DistributionChangeCheck(&out, 6, 0);
 277   defocusLikelihoodInit(&likelihood,&out,1);
 278
 279   DefocusEstimationLoop(info,&out,&ctfinfo,&ctf,&intensity,&likelihood,6,1);
 280   DEBUGPRINT("defocus refinement finish\n");
 281   ctfforbayescheckIntensityOfMaxPosterior(ctf,out.posterior,out.numDistribution,&ctfinfo,intensity.MTFmode,intensity.Ainmode,1,1);
 282
 283   eosBayesPosteriorWrite(&out,0);
 284
 285   ctfforbayesFree(&ctf,&intensity,&likelihood,out.numLikelihood);
 286
 287   eosBayesFree(&out);
 288
 289   exit(EXIT_SUCCESS);
 290 }
 291
 292 void DefocusEstimationLoop(ctfDeterminationByBayesInfo info, eosBayes* out, ctfInfo* ctfinfo, ctfforbayes* ctf, intensityforbayes* intensity, defocuslikelihoodforbayes* likelihood, int count,int mode){
 293   int i;
 294
 295   defocusLikelihoodInit(likelihood,out,mode);
 296   DEBUGPRINT("In DefocusEstimationLoop\n");
 297   for(i=0;i<info.flagIn4;i++){
 298
 299         intensityRead(info.In4[i],intensity,i,0);
 300         DEBUGPRINT("ctf start\n");
 301         ctfFunctionforbayes(ctf,intensity,out,ctfinfo,likelihood,count,1);
 302
 303         DEBUGPRINT("estimate start\n");
 304         eosBayesEstimation(out,0);
 305         DistributionChangeCheck(out,count,1);
 306
 307         DEBUGPRINT("change start\n");
 308         eosBayesPosteriortoPrior(out,0);
 309
 310   }
 311
 312 }
 313
 314 void
 315 additionalUsage()
 316 {
 317   fprintf(stderr, "----- Additional Usage -----\n");
 318 }
 319