src/Objects/DataManip/ctfInfo/src/org/lctfDetermination.c.12_6

   1 /*
   2 # %M% %Y% %I%
   3 # The latest update : %G% at %U%
   4 #
   5 #%Z% lctfDetermination ver %I%
   6 #%Z% Created by
   7 #%Z%
   8 #%Z% Usage : lctfDetermination
   9 #%Z% Attention
  10 #%Z%
  11 */
  12 static char __sccs_id[] = "%Z%lctfDetermination ver%I%; Date:%D% %Z%";
  13 #include <stdio.h>
  14 #include <stdlib.h>
  15 #include <math.h>
  16 #include "./lctfDetermination.h"
  17 #define DEBUG
  18 #include "genUtil.h"
  19 #include "Vector.h"
  20 #include "nr2.h"
  21
  22 /*
  23     Henderson-like Min-Max Method
  24                 0.9<CTF*CTF: +Scattering : Score
  25                 CTF*CTF<0.1: -Scattering : Penalty
  26 */
  27
  28 void
  29 lctfDeterminationbyMinMaxMethods(ctfInfo* res, mrcImage* mrc, ctfInfo* ini, long mode)
  30 {
  31     floatVector* scatter;
  32     floatVector* spacing;
  33     floatVector* ctf;
  34     long         i;
  35     float        E;
  36         float        maxE;
  37     ctfInfo      tmp;
  38
  39     scatter = lmrcFSInfoScatteringAngularDistributionAverage(mrc);
  40     spacing = lmrcFSInfoSpacing(mrc);
  41         DEBUGPRINT2("%f %f\n", scatter->data[1], spacing->data[1]);
  42     ctf     = floatVectorInit(NULL, spacing->size);
  43     tmp     = *ini;
  44     *res    = *ini;
  45
  46     maxE = 0;
  47
  48     for(tmp.defocus=0; tmp.defocus<2*ini->defocus; tmp.defocus+=ini->defocus/100.0) {
  49         for(i=0; i<spacing->size; i++) {
  50             ctf->data[i] = ctfFunction(&tmp, spacing->data[i], mode);
  51         }
  52                 DEBUGPRINT1("%f\n", ctf->data[1]);
  53         E = 0.0;
  54         for(i=0; i<spacing->size; i++) {
  55             if(tmp.CutLow < spacing->data[i]
  56              &&spacing->data[i] < tmp.CutHigh ) {
  57                 if(SQR(ctf->data[i])>0.9) {
  58                     E += scatter->data[i];
  59                 } else if(SQR(ctf->data[i])<0.1) {
  60                     E -= scatter->data[i];
  61                 }
  62             }
  63         }
  64         if(maxE<E) {
  65                         fprintf(stdout, "Max: %f %f %f\n", tmp.defocus, E, maxE);
  66                         maxE = E;
  67             res->defocus = tmp.defocus;
  68         }
  69         fprintf(stdout, "%10.6e %10.6e\n", tmp.defocus, E);
  70     }
  71 }
  72
  73 /*
  74         Fitting Method:
  75
  76     F(R) = S(R)*S(R) + N(R)*N(R) + N2(R)*N2(R)
  77         , where
  78                 S(R)  = I0*|CTF(R)|*Env(R)*ME(R)*Vib(R)*MTF(R) ,
  79                 N(R)  = I0*White*MTF(R) and
  80                 N2(R) = I0*White2
  81
  82         Fitted to F(R)
  83         A   : Normalization Factor
  84         CTF :   -sin(Kai)-a*cos(Kai)
  85                 Kai = M_PI*dF*lambda*R*R - 0.5*M_PI*Cs*lambda*lambda*lambda*R*R*R*R
  86                 ENV(Cs) : exp(-SQR(M_PI*Cs*lambda*lambda*R*R*R-M_PI*dF*R)*SQR(Ain)/log(2))
  87                 ENV(Cc) : exp(-SQR(M_PI*lambda*R*R*Cc*dE/E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2))))
  88         Me  : Table
  89         SN  : constant (White Noise)
  90         B   : constant (MTF) SingleExp/Lorenzian/Both
  91
  92     a[1]  : I0  : normarizing factor
  93     a[2]  : dF : defocus (underfocus is positive)
  94     a[3]  : Cs : Cs
  95     a[4]  : a  : the ratio of amplitude contrast to phase contrast
  96     a[5]  : Ai : Illumination angle
  97     a[6]  : Cc : Cc
  98     a[7]  : dE : dE
  99     a[8]  : V  : Vibration : exp(-V*V*R*R/2)
 100     a[9]  : White  : signal to noise ratio
 101     a[10] : White2 : signal to noise ratio
 102     a[11] : B  : Dumping Factor of MTF: exp(-B*R)
 103     a[12] : E  : Energy : fixed
 104     a[13] : E0 : Energy : fixed
 105     a[14] : lambda  : Wave Length : fixed
 106     a[15] : Me : Scattering Factor : fixed
 107
 108 */
 109 static ctfInfo __lctfDetermineCTF;
 110
 111
 112 void
 113 angularDistributionFunctionOfSignalAndNoise(float x, float p[], float* y, float dyda[], int na)
 114 {
 115     float R, I0, dF, Cs, a, Ai, Cc, dE, V, White, White2, B,  E, E0, lambda, Me;
 116         float chi, dchidR, dchi, spread, dspread, CTF, dCTF, Env1, Env2, Env, Vib, MTF;
 117         float S;
 118         float N;
 119         float N2;
 120         float F;
 121         float Mag;
 122
 123         Mag = p[16];
 124     R  = x*Mag;
 125     I0 = p[1];
 126     dF = p[2];
 127     Cs = p[3];
 128     a  = p[4];
 129     Ai = p[5];
 130     Cc = p[6];
 131     dE = p[7];
 132     V  = p[8];
 133     White  = p[9];
 134     White2 = p[10];
 135     B  = p[11];
 136         /* */
 137         if(__lctfDetermineCTF.flagMolecEnvTable) {
 138                 Me = lmolecularEnvelopeDataGet(&(__lctfDetermineCTF.MolecEnvTable), R*1e-10, 0);
 139         } else {
 140                 Me = exp(-p[12]*p[12]*R*R/2.0);
 141         }
 142         /* fixed */
 143     E  = p[13];
 144     E0 = p[14];
 145     lambda = p[15];
 146
 147         /* Wave Length */
 148     chi    = M_PI*(dF*lambda - 0.5*Cs*lambda*lambda*lambda*R*R)*R*R;
 149     dchidR = 2*M_PI*(dF*lambda - Cs*lambda*lambda*lambda*R*R)*R;
 150         dchi   = M_PI*(Cs*lambda*lambda*R*R - dF)*R;
 151         spread = M_PI*lambda*R*R*Cc*dE/E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2));
 152     CTF = -sin(chi) - a*cos(chi);
 153         Env1 = exp(-SQR(dchi)*SQR(Ai)/log(2));
 154         Env2 = exp(-SQR(spread));
 155         Env  = Env1*Env2;
 156         Vib  = exp(-V*V*R*R/2.0);
 157         MTF  = exp(-B*R);
 158
 159         S  = I0*Me*CTF*Env*Vib*MTF;
 160         N  = I0*White*MTF;
 161         N2 = I0*White2;
 162
 163     F  = SQR(S) + SQR(N) + SQR(N2);
 164         *y = F;
 165
 166
 167     dyda[1] = 2*F/I0;
 168     dyda[2] = 2*S*I0*Me*(-cos(chi)+a*sin(chi))*(M_PI*lambda*R*R)*Env*Vib*MTF
 169                  +2*S*S*(-2)*dchi*(-M_PI*R)*SQR(Ai)/log(2);
 170     dyda[3] = 2*S*I0*Me*(-cos(chi)+a*sin(chi))*(M_PI*(-0.5)*lambda*lambda*lambda*R*R*R*R)*Env*Vib*MTF
 171                  +2*S*S*(-2)*dchi*(M_PI*lambda*lambda*R*R*R)*SQR(Ai)/log(2);
 172     dyda[4] = 2*S*I0*Me*(-cos(chi))*Env*Vib*MTF;
 173     dyda[5] = 2*S*S*(-2)*SQR(dchi)*Ai/log(2);
 174     dyda[6] = 2*S*S*(-2)*spread*M_PI*lambda*R*R*   dE/E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2));
 175     dyda[7] = 2*S*S*(-2)*spread*M_PI*lambda*R*R*Cc   /E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2));
 176     dyda[8] =   S*S*(-1)*V*R*R;
 177     dyda[9] = 2*N *N /White;
 178     dyda[10]= 2*N2*N2/White2;
 179     dyda[11]= 2*(S*S+N*N)*(-R);
 180         /* fixed */
 181         if(__lctfDetermineCTF.flagMolecEnvTable) {
 182         dyda[12] = 0;
 183         } else {
 184                 dyda[12] = 2*S*S*(-p[12]*R*R);
 185         }
 186         DEBUGPRINT6("%e %e %e %e %e %e ", dyda[1], dyda[2], dyda[3], dyda[4], dyda[5], dyda[6]);
 187         DEBUGPRINT6("%e %e %e %e %e %e \n", dyda[7], dyda[8], dyda[9], dyda[10], dyda[11], dyda[12]);
 188     dyda[13]= 0.0;
 189     dyda[14]= 0.0;
 190     dyda[15]= 0.0;
 191
 192         if(__lctfDetermineCTF.flagMolecEnvTable) {
 193                 dyda[16]= x*
 194                                   (2*S*
 195                                   (I0*Me*((-cos(chi)+a*sin(chi))*dchidR)*Env*Vib*MTF
 196                                  + I0*Me*CTF*((-2)*dchi*SQR(Ai)/log(2)*M_PI*(3*Cs*lambda*lambda*R - dF)*Env1)*Vib*MTF
 197                                  + I0*Me*CTF*((-2)*spread*2*M_PI*lambda*R*Cc*dE/E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2))*Env2)*Vib*MTF
 198                                  + I0*Me*CTF*Env*(-V*V*R*Vib)*MTF
 199                                  + I0*Me*CTF*Env*Vib*(-B*MTF))
 200                                  +(2*N*(-B*N)));
 201         } else {
 202                 dyda[16]= x*
 203                                   (2*S*
 204                                   (I0*(-V*V*R*Me)*CTF*Env*Vib*MTF
 205                                  + I0*Me*((-cos(chi)+a*sin(chi))*dchidR)*Env*Vib*MTF
 206                                  + I0*Me*CTF*((-2)*dchi*SQR(Ai)/log(2)*M_PI*(3*Cs*lambda*lambda*R*R - dF)*Env1)*Env2*Vib*MTF
 207                                  + I0*Me*CTF*Env1*((-2)*spread*2*M_PI*lambda*R*Cc*dE/E*((1+E/E0)/(1+E/E0/2))/4/sqrt(log(2))*Env2)*Vib*MTF
 208                                  + I0*Me*CTF*Env*(-V*V*R*Vib)*MTF
 209                                  + I0*Me*CTF*Env*Vib*(-B*MTF))
 210                                  +(2*N*(-B*N)));
 211         }
 212 }
 213
 214 void
 215 lctfDeterminationbyFittingMethods(ctfInfo* res, ctfInfo* var, mrcImage* mrc, ctfInfo* ini, long mode)
 216 {
 217
 218     floatVector* scatter;
 219     floatVector* spacing;
 220     floatVector*  ctf;
 221     floatVector*  sig;
 222     long i, j, imin, imax, n, flag;
 223     float E, maxE;
 224     ctfInfo tmp;
 225         int count;
 226     float* a;
 227     int* ia;
 228     int ma;
 229     float** covar;
 230     float** alpha;
 231     float chisq;
 232     float alambda, oldalambda;
 233
 234     scatter = lmrcFSInfoScatteringAngularDistributionAverage(mrc);
 235     sig     = lmrcFSInfoScatteringAngularDistributionSD(mrc);
 236     spacing = lmrcFSInfoSpacing(mrc);
 237     ctf     = floatVectorInit(NULL, spacing->size);
 238     imin = 0; imax = spacing->size - 1; n = 0;
 239     for(i=0; i<spacing->size; i++) {
 240         if(ini->CutLow  <= spacing->data[i]
 241          &&ini->CutHigh >= spacing->data[i]) {
 242             if(0==n) {
 243                 imin = i;
 244             }
 245             n++;
 246         }
 247         spacing->data[i] = spacing->data[i]*1e10 ;
 248         DEBUGPRINT4("%d: %g %g %g\n", i, spacing->data[i], scatter->data[i], sig->data[i]);
 249     }
 250         if(n<3) {
 251                 fprintf(stderr, "mrcImage is too small : %d!!\n", n);
 252                 exit(EXIT_FAILURE);
 253         }
 254     DEBUGPRINT4("Low: %g High: %g ::: n: %d, imin: %d\n", ini->CutLow, ini->CutHigh, n, imin);
 255     tmp  = *ini;
 256     *res = *ini;
 257
 258     ma = 16;
 259     a  = vector(1, ma);
 260     ia = ivector(1, ma);
 261     covar = matrix(1, ma, 1, ma);
 262     alpha = matrix(1, ma, 1, ma);
 263
 264         for(i=1; i<=ma; i++) {
 265                 a[i]  = 0;
 266                 ia[i] = 0;
 267                 for(j=1; j<=ma; j++) {
 268                         covar[i][j] = 0;
 269                         alpha[i][j] = 0;
 270                 }
 271         }
 272         /* Parameter */
 273     a[1] = ini->I0;                 ia[1] = 1;
 274     a[2] = ini->defocus*1e-10;      ia[2] = 1;
 275     a[3] = ini->Cs*1e-3;            ia[3] = 1;
 276     a[4] = ini->ratioOfAmpToPhase;  ia[4] = 1;
 277     a[5] = ini->Ain*1e-3;           ia[5] = 1;
 278     a[6] = ini->Cc*1e-3;            ia[6] = 1;
 279     a[7] = ini->dE;                 ia[7] = 1;
 280     a[8] = ini->BofVibration*1e-10; ia[8] = 1;
 281     a[9] = ini->WhiteNoise;         ia[9] = 1;
 282     a[10]= ini->WhiteNoise2;        ia[10] = 1;
 283     a[11]= ini->BofMTF*1e-10;       ia[11] = 1;
 284         if(ini->flagMolecEnvTable) {
 285         a[12]= 0;   ia[12] = 0;
 286         } else {
 287         a[12]= ini->MolecEnv*1e-10; ia[12] = 1;
 288         }
 289         /* fixed */
 290     a[13] = ini->kV*1e3;                 ia[13] = 0;
 291     a[14] = 511*1e3;                     ia[14] = 0;
 292     a[15] = wavelengthOfElectron(ini->kV*1e3); ia[15] = 0;
 293         fprintf(stderr, "lambda: %g\n", a[15]);
 294         /* Parameter */
 295         a[16] = ini->Magnification;     ia[16] = 1;
 296
 297         __lctfDetermineCTF = *ini;
 298
 299     alambda = -1;
 300         oldalambda = 1e-6;
 301         for(count=0; count<20; count++) {
 302                 fprintf(stdout, "count: %d\n", count); fflush(stdout);
 303         mrqmin(spacing->data+imin-1, scatter->data+imin-1, sig->data+imin-1, n, a, ia, ma, covar, alpha, &chisq, angularDistributionFunctionOfSignalAndNoise, &alambda);
 304                 fprintf(stdout, "chisq: %15.6g\n", chisq);
 305                 fprintf(stdout, "alambda: %15.6g\n", alambda);
 306                 for(i=1; i<=ma; i++) {
 307                         fprintf(stdout, "a[%d]: %15.6g\n", i, a[i]);
 308                 }
 309                 if(alambda<1e-6) {
 310                         break;
 311                 } else if(oldalambda<alambda) {
 312                         count--;
 313                 }
 314                 oldalambda = alambda;
 315                 fflush(stdout);
 316         }
 317         fprintf(stdout, "Fitting End: %15.6g\n", chisq);
 318     res->I0                = a[1];
 319     res->defocus           = a[2]*1e10;
 320     res->Cs                = a[3]*1e+3;
 321     res->ratioOfAmpToPhase = a[4];
 322     res->Ain               = a[5]*1e+3;
 323     res->Cc                = a[6]*1e+3;
 324     res->dE                = a[7];
 325     res->BofVibration      = a[8]*1e10;
 326     res->WhiteNoise        = a[9];
 327     res->WhiteNoise2       = a[10];
 328     res->BofMTF            = a[11]*1e10;
 329         if(ini->flagMolecEnvTable) {
 330         } else {
 331         res->MolecEnv      = a[12]*1e10;
 332         }
 333         res->Magnification     = a[16];
 334
 335         alambda = 0;
 336     mrqmin(spacing->data+imin-1, scatter->data+imin-1, sig->data+imin-1, n, a, ia, ma, covar, alpha, &chisq, angularDistributionFunctionOfSignalAndNoise, &alambda);
 337         fprintf(stdout, "chisq: %15.6g\n", chisq);
 338         fprintf(stdout, "alambda: %15.6g\n", alambda);
 339         for(i=1; i<=ma; i++) {
 340                 fprintf(stdout, "a[%d]: %15.6g\n", i, a[i]);
 341         }
 342     var->I0                = sqrt(covar[1][1]);
 343     var->defocus           = sqrt(covar[2][2])*1e10;
 344     var->Cs                = sqrt(covar[3][3])*1e+3;
 345     var->ratioOfAmpToPhase = sqrt(covar[4][4]);
 346     var->Ain               = sqrt(covar[5][5])*1e+3;
 347     var->Cc                = sqrt(covar[6][6])*1e+3;
 348     var->dE                = sqrt(covar[7][7]);
 349     var->BofVibration      = sqrt(covar[8][8])*1e10;
 350     var->WhiteNoise        = sqrt(covar[9][9]);
 351     var->WhiteNoise2       = sqrt(covar[10][10]);
 352     var->BofMTF            = sqrt(covar[11][11])*1e10;
 353         if(ini->flagMolecEnvTable) {
 354         } else {
 355         res->MolecEnv      = sqrt(covar[12][12])*1e10;
 356         }
 357         var->Magnification     = sqrt(covar[16][16]);
 358 }
 359
 360 /*
 361         Differential Method
 362                 ThonRing Point: ScatteringDifferential  :0      score+
 363                                                         :+or-   score-
 364 */
 365 /*
 366         LimitValue
 367                 determination      MaximumValue         :return         +1
 368                                 or minimumValue         :return         -1
 369                                    NotLimitValue        :return          0
 370 */
 371
 372 static int  __LimitValue(floatVector* function,long i,long around)
 373 {
 374         float   before = 0.0,after = 0.0;
 375         long    j;
 376         int     flag;
 377
 378         j = i-around;
 379         //while(before == 0.0){
 380                 before  =function->data[i]-function->data[j];
 381         //      j -= 1;
 382         //}
 383
 384         j = i+around;
 385         //while(after == 0.0){
 386                 after   =function->data[j]-function->data[i];
 387         //      j += 1;
 388         //}
 389
 390         if(before > 0.0){
 391                 if(after <= 0){
 392                         flag = 1;
 393                 }
 394                 else{
 395                         flag = 0;
 396                 }
 397         }
 398         else if(before < 0.0){
 399                 if(after >= 0){
 400                         flag = -1;
 401                 }
 402                 else{
 403                         flag = 0;
 404                 }
 405         }
 406         else{
 407                 if(after > 0.0){
 408                         flag = -1;
 409                 }
 410                 else{
 411                         flag = +1;
 412                 }
 413         }
 414         return(flag);
 415 }
 416
 417 static int  __SquareLimitValue(floatVector* function,long i,long around)
 418 {
 419         float   before = 0.0,after = 0.0;
 420         long    j;
 421         int     flag;
 422
 423         j = i-around;
 424         //while(before == 0.0){
 425                 before  = SQR(function->data[i])-SQR(function->data[j]);
 426         //      j -= 1;
 427         //}
 428
 429         j = i+around;
 430         //while(after == 0.0){
 431                 after   = SQR(function->data[j])-SQR(function->data[i]);
 432         //      j += 1;
 433         //}
 434         //fprintf(stdout,"%ld %f %f ",i,before,after);
 435
 436
 437         if(before > 0.0){
 438                 if(after <= 0){
 439                         flag = 1;
 440                 }
 441                 else{
 442                         flag = 0;
 443                 }
 444         }
 445         else if(before < 0.0){
 446                 if(after >= 0){
 447                         flag = -1;
 448                 }
 449                 else{
 450                         flag = 0;
 451                 }
 452         }
 453         else{
 454
 455                 if(after > 0.0){
 456                         flag = -1;
 457                 }
 458                 else{
 459                         flag = +1;
 460                 }
 461         }
 462         //fprintf(stdout,"%ld\n",flag);
 463         return(flag);
 464 }
 465
 466
 467 static float __Median(floatVector* function, long i, long around){
 468
 469         float sequence[around*2 + 1];
 470         float tmp;
 471         long  j,k;
 472
 473         k = 0;
 474         for(j=i-around;j<=i+around;j++){
 475                 sequence[k] = function->data[j];
 476                 k++;
 477         }
 478
 479         for(j=0;j<around*2+1;j++){
 480                 for(k=j+1;k<around*2+1;k++){
 481                         if(sequence[j]>sequence[k]){
 482                                 tmp = sequence[j];
 483                                 sequence[j] = sequence[k];
 484                                 sequence[k] = tmp;
 485                         }
 486                 }
 487         }
 488         return(sequence[around]);
 489 }
 490
 491 static float __Average(floatVector* function,long i,long around){
 492
 493         long  j;
 494         float sum = 0.0;
 495
 496         for(j=i-around;j<=i+around;j++){
 497                 sum += function->data[j];
 498         }
 499
 500         return(sum/(around*2+1));
 501 }
 502
 503
 504 static long __LongAbs(long x){
 505
 506         if(x <=  0){
 507                 return(-x);
 508         }
 509         else{
 510                 return(x);
 511         }
 512 }
 513 static float __AverageWithWeight(floatVector* function,long i,long around){
 514
 515         long  j;
 516         float sum = 0.0;
 517         float n = 0.0;
 518
 519         for(j=i-around;j<=i+around;j++){
 520                 sum += function->data[j]/(__LongAbs(i-j)+1);
 521                 n   += (float)1/( __LongAbs(i-j) + 1);
 522         }
 523         return( sum/n );
 524
 525 }
 526
 527
 528 /*Smoothing is CurveSmoothing
 529         mode 1: for MedianMethod
 530         mode 2: for AverageMethod
 531         mode 3: for AverageMethod with Weight
 532 */
 533 void __Smoothing(floatVector* function,floatVector* SmoothCurve ,long around,long mode){
 534
 535         long i;
 536         long aroundtmp;
 537
 538         aroundtmp = around;
 539
 540         for(i=0;i<function->size;i++){
 541                 aroundtmp = around;
 542                 if(i!=0 && i!=function->size-1){
 543                         if(i<around){
 544                                 aroundtmp = i;
 545                         }
 546                         else if(i>function->size-1-around){
 547                                 aroundtmp = function->size-1-i;
 548                         }
 549
 550                         switch(mode){
 551                                 case 1:{
 552                                         SmoothCurve->data[i] = __Median(function,i,aroundtmp);
 553                                         break;
 554                                 }
 555                                 case 2:{
 556                                         SmoothCurve->data[i] = __Average(function,i,aroundtmp);
 557                                         break;
 558                                 }
 559                                 case 3:{
 560                                         SmoothCurve->data[i] = __AverageWithWeight(function,i,aroundtmp);
 561                                         break;
 562                                 }
 563                         }
 564                 }
 565
 566                 else{
 567                         SmoothCurve->data[i] = function->data[i];
 568                 }
 569         }
 570
 571 }
 572
 573
 574 void
 575 lctfDeterminationbyDifferentialMethods(ctfInfo* res, mrcImage* mrc, ctfInfo* ini,long* MaxThonRing, long mode)
 576 {
 577         floatVector* scatter;
 578         floatVector* spacing;
 579         floatVector* ctf;
 580         floatVector* MedianCurve;
 581         long         i;
 582         long         Differentialaround;
 583         long         Smootharound;
 584         long         Smoothmode;
 585         long         ThonRing;
 586         float        E;
 587         float        maxE=0;
 588         float        ME;                        /*M......with MedianCurve*/
 589         float        MmaxE=0;
 590         ctfInfo      tmp;
 591         ctfInfo      Mres;
 592         //float        Mres;
 593
 594         Differentialaround = 1;
 595         Smootharound       = 2;
 596         Smoothmode         = 1;
 597
 598         scatter = lmrcFSInfoScatteringAngularDistributionAverage(mrc);
 599         spacing = lmrcFSInfoSpacing(mrc);
 600         DEBUGPRINT2("%f %f\n", scatter->data[1], spacing->data[1]);
 601         ctf     = floatVectorInit(NULL, spacing->size);
 602         tmp     = *ini;
 603         *res    = *ini;
 604         Mres   = *ini;
 605
 606         fprintf(stdout,"\nDifferentialMethod\n");
 607         MedianCurve = floatVectorInit(NULL,spacing->size);
 608         for(tmp.defocus = 0.0; tmp.defocus<2*ini->defocus; tmp.defocus+=ini->defocus/1000.0) {
 609                 for(i=0;i<spacing->size;i++){
 610                         ctf->data[i] = ctfFunction(&tmp, spacing->data[i],mode);
 611                 }
 612                 __Smoothing(scatter,MedianCurve,Smootharound,Smoothmode);
 613                 DEBUGPRINT1("%f\n", ctf->data[1]);
 614                 E = 0.0;
 615                 ME= 0.0;
 616                 ThonRing = 0;
 617                 for(i=(0+Differentialaround); i<(spacing->size-Differentialaround); i++) {
 618                         if(tmp.CutLow < spacing->data[i]
 619                         &&spacing->data[i] < tmp.CutHigh ) {
 620                                 /*DifferentialMethod*/
 621                                 if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
 622                                         ThonRing += 1;
 623                                         if(__LimitValue(scatter,i,Differentialaround) == -1){
 624                                                 E+=1.0;
 625                                         }
 626                                         else{
 627                                                 E-=1.0;
 628                                         }
 629
 630                                         if(__LimitValue(MedianCurve,i,Differentialaround) == -1){
 631                                                 ME+=1.0;
 632                                         }
 633                                         else{
 634                                                 ME-=1.0;
 635                                         }
 636
 637                                 }
 638                         }
 639                 }
 640                 if(maxE<E) {
 641                         fprintf(stdout, "Max: %10.6e %f %f\n", tmp.defocus, E, maxE);
 642                         maxE = E;
 643                         res->defocus = tmp.defocus;
 644                 }
 645
 646                 if(MmaxE<ME){
 647                         //fprintf(stdout, "MMAX: %f %f %f\n", tmp.defocus,ME,MmaxE);
 648                         MmaxE = ME;
 649                         Mres.defocus         = tmp.defocus;
 650                         *MaxThonRing = ThonRing;
 651                 }
 652                 //fprintf(stdout, "%10.6e %3ld  %3.0f %3.0f\n", tmp.defocus,ThonRing, E,ME);
 653         }
 654         fprintf(stdout,"defocus(not median)= %10.6f (median)= %10.6f\n",res->defocus,Mres.defocus);
 655
 656         res->defocus = Mres.defocus;
 657         /*for smooth test*/
 658         /*for(i=0;i<spacing->size;i++){
 659                 ctf->data[i]  = ctfFunction(&Mres, spacing->data[i],mode);
 660         }
 661         fprintf(stdout,"\n\ndefocus:%f\n",Mres.defocus);
 662         fprintf(stdout,"    spacing:    SQR(ctf):     scatter:\n");
 663         for(i=0;i<spacing->size;i++){
 664
 665                 if(tmp.CutLow < spacing->data[i]
 666                 &&spacing->data[i] < tmp.CutHigh ) {
 667                         fprintf(stdout,"%10.6e %10.6e %10.6e ",spacing->data[i],SQR(ctf->data[i]),scatter->data[i]);
 668                         if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
 669                                 fprintf(stdout,"ring\n");
 670                         }
 671                         else{
 672                                 fprintf(stdout,"\n");
 673                         }
 674                 }
 675         }*/
 676 }
 677
 678 void __MinSquare(floatVector* functionX,floatVector* functionY,floatVector* MinSquare,long around){
 679         long i;
 680         long j;
 681         long tmpAround;
 682         float x;
 683         float y;
 684         float xy;
 685         float xx;
 686         long n;
 687
 688         tmpAround = around;
 689         n = around*2+1;
 690         for(i=around;i<functionX->size-around;i++){
 691                 tmpAround = around;
 692
 693                 if(i<around){
 694                         tmpAround = i;
 695                 }
 696                 else if(functionX->size-1-i<around){
 697                         tmpAround = functionX->size-1-i;
 698                 }
 699
 700                 if(i!=0 || i!=functionX->size-1){
 701                         x = 0.0;
 702                         y = 0.0;
 703                         xx = 0.0;
 704                         xy = 0.0;
 705                         for(j=i-around;j<=i+around;j++){
 706                                 x += functionX->data[j];
 707                                 y += functionY->data[j];
 708                                 xy += functionX->data[j]*functionY->data[j];
 709                                 xx += functionX->data[j]*functionX->data[j];
 710                         }
 711                         MinSquare->data[i] = (n*xy-x*y)/(n*xx-x*x);
 712                 }
 713                 else{
 714                         MinSquare->data[i] = 0.0;
 715                 }
 716         }
 717 }
 718
 719 void __MinimumByMinSquare(floatVector* functionX,floatVector* functionY,floatVector* Minimum,long around)
 720 {
 721         floatVector* MinSquare;
 722         long i;
 723         long j;
 724         long k;
 725
 726         MinSquare = floatVectorInit(NULL,functionX->size);
 727         __MinSquare(functionX,functionY,MinSquare,around);
 728
 729         for(i=0;i<functionX->size;i++){
 730                 if(MinSquare->data[i] > 0.0){
 731                         MinSquare->data[i] = 1.0;
 732                 }
 733                 else if(MinSquare->data[i] < 0.0){
 734                         MinSquare->data[i] = -1.0;
 735                 }
 736
 737         }
 738
 739         for(i=1;i<functionX->size;i++){
 740                 switch ((int)MinSquare -> data[i-1]){
 741                         case (int)-1:{
 742                                 switch ((int)MinSquare->data[i]){
 743                                         case (int)1:{
 744                                                 Minimum->data[i] = 1.0;
 745                                                 break;
 746                                         }
 747                                         case (int)-1:{
 748                                                 Minimum->data[i] = 0.0;
 749                                                 break;
 750                                         }
 751                                         case (int)0:{
 752                                                 j = 1;
 753                                                 while(MinSquare->data[i+j] == 0.0 ){
 754                                                         if(i+j == functionX->size){
 755                                                                 break;
 756                                                         }
 757                                                         j++;
 758                                                 }
 759
 760                                                 if(MinSquare->data[i+j] > 0){
 761                                                         Minimum->data[i] == 1.0;
 762                                                 }
 763                                                 else if(MinSquare->data[i+j] < 0){
 764                                                         Minimum->data[i] == 0.0;
 765                                                 }
 766                                                 break;
 767                                         }
 768                                 }
 769                         }
 770
 771
 772                         case (int)1:{
 773                                 Minimum->data[i] == 0.0;
 774                                 break;
 775                         }
 776
 777                         case (int)0:{
 778                                 j = -1;
 779                                 while(MinSquare->data[i+j] == 0.0){
 780                                         if(i+j == 0){
 781                                                 break;
 782                                         }
 783                                         j--;
 784                                 }
 785
 786                                 if(MinSquare->data[i+j] > 0.0){
 787                                         Minimum->data[i] == 0.0;
 788                                 }
 789                                 else if(MinSquare->data[i+j] <= 0.0){
 790                                         k = 1;
 791                                         while(MinSquare->data[i+k] == 0.0){
 792                                                 if(i+k == functionX->size){
 793                                                         break;
 794                                                 }
 795                                                 k++;
 796                                         }
 797
 798                                         if(MinSquare->data[i+k] >= 0.0){
 799                                                 Minimum->data[i] = 1.0;
 800                                         }
 801                                         else{
 802                                                 Minimum->data[i] = 0.0;
 803                                         }
 804                                 }
 805                                 break;
 806                         }
 807                 }
 808         }
 809 }
 810
 811 void
 812 lctfDeterminationbyMinSquareMethods(ctfInfo* res, mrcImage* mrc, ctfInfo* ini,long* MaxThonRing, long mode)
 813 {
 814         floatVector* scatter;
 815         floatVector* spacing;
 816         floatVector* ctf;
 817         floatVector* ctfSQR;
 818         floatVector* MedianCurve3;
 819         floatVector* MedianCurve5;
 820         floatVector* MedianCurve7;
 821         floatVector* MinSquare;
 822         floatVector* Minimum;
 823         floatVector* ThonRing;
 824         floatVector* ctfMinSquare;
 825         long         i;
 826         long         Differentialaround;
 827         long         Smootharound;
 828         long         Smoothmode;
 829         long         ThonRingSize;
 830         float        E;
 831         float        maxE=0;
 832         float        ME;                        /*M......with MedianCurve*/
 833         float        MmaxE=0;
 834         ctfInfo      tmp;
 835         ctfInfo      Mres;
 836
 837         Differentialaround = 1;
 838         Smoothmode         = 1;
 839
 840         scatter = lmrcFSInfoScatteringAngularDistributionAverage(mrc);
 841         spacing = lmrcFSInfoSpacing(mrc);
 842         DEBUGPRINT2("%f %f\n", scatter->data[1], spacing->data[1]);
 843         ctf     = floatVectorInit(NULL, spacing->size);
 844         ctfSQR  = floatVectorInit(NULL, spacing->size);
 845         ThonRing= floatVectorInit(NULL, spacing->size);
 846         ctfMinSquare = floatVectorInit(NULL,spacing->size);
 847         tmp     = *ini;
 848         *res    = *ini;
 849         Mres   = *ini;
 850
 851         fprintf(stdout,"\nDifferentialMethod\n");
 852         MedianCurve3 = floatVectorInit(NULL,spacing->size);
 853         MedianCurve5 = floatVectorInit(NULL,spacing->size);
 854         MedianCurve7 = floatVectorInit(NULL,spacing->size);
 855         MinSquare    = floatVectorInit(NULL,spacing->size);
 856         Minimum      = floatVectorInit(NULL,spacing->size);
 857
 858         __Smoothing(scatter,MedianCurve3,1,1);
 859         __Smoothing(scatter,MedianCurve5,2,1);
 860         __Smoothing(scatter,MedianCurve7,3,1);
 861
 862         __MinimumByMinSquare(spacing,MedianCurve5,Minimum,2);
 863         __MinSquare(spacing,MedianCurve5,MinSquare,2);
 864
 865         for(tmp.defocus = 0.0; tmp.defocus<2*ini->defocus; tmp.defocus+=ini->defocus/1000.0) {
 866                 for(i=0;i<spacing->size;i++){
 867                         ctf->data[i] = ctfFunction(&tmp, spacing->data[i],mode);
 868                 }
 869                 DEBUGPRINT1("%f\n", ctf->data[1]);
 870                 E = 0.0;
 871                 ME= 0.0;
 872                 ThonRingSize = 0;
 873                 for(i=(0+Differentialaround); i<(spacing->size-Differentialaround); i++) {
 874                         if(tmp.CutLow < spacing->data[i]
 875                         &&spacing->data[i] < tmp.CutHigh ) {
 876
 877                                 /*DifferentialMethod*/
 878                                 if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
 879                                         ThonRingSize += 1;
 880                                         if(/*__LimitValue(MedianCurve5,i,Differentialaround)*/Minimum->data[i] == 1){
 881                                                 ME+=1.0;
 882                                         }
 883                                         else{
 884                                                 ME+=-1.0;
 885                                         }
 886
 887                                 }
 888
 889
 890                         }
 891                 }
 892
 893                 if(MmaxE<ME){
 894                         fprintf(stdout, "MMAX: %f %f %f\n", tmp.defocus,ME,MmaxE);
 895                         MmaxE = ME;
 896                         res->defocus         = tmp.defocus;
 897                         *MaxThonRing = ThonRingSize;
 898                 }
 899                 //fprintf(stdout, "%10.6e %3ld  %3.0f %3.0f\n", tmp.defocus,ThonRingSize, E,ME);
 900         }
 901         fprintf(stdout,"defocus(median)= %10.6f\n",res->defocus);
 902
 903         /*for smooth test*/
 904         for(i=0;i<spacing->size;i++){
 905                 ctf->data[i]  = ctfFunction(res, spacing->data[i],mode);
 906                 ctfSQR->data[i] = SQR(ctf->data[i]);
 907         }
 908         __MinimumByMinSquare(spacing,ctfSQR,ThonRing,1);
 909         __MinSquare(spacing,ctfSQR,ctfMinSquare,1);
 910         fprintf(stdout,"\n\ndefocus:%f\n",res->defocus);
 911         fprintf(stdout,"    spacing:    SQR(ctf):     scatter:     Median:  MinSquare:    Minimum:\n");
 912         for(i=0;i<spacing->size;i++){
 913
 914                 if(tmp.CutLow < spacing->data[i]
 915                 &&spacing->data[i] < tmp.CutHigh ) {
 916                         fprintf(stdout,"%10.6e %10.6e %10.6e %10.6e %+10.6e %f %+10.6e %f ",spacing->data[i],SQR(ctf->data[i]),scatter->data[i],MedianCurve5->data[i],MinSquare->data[i],Minimum->data[i],ctfMinSquare->data[i],ThonRing->data[i]);
 917                         if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
 918                                 fprintf(stdout,"ring\n");
 919                         }
 920                         else{
 921                                 fprintf(stdout,"\n");
 922                         }
 923                 }
 924         }
 925 }
 926
 927
 928 /*ctfDeterminationbyCovarianceMethod*/
 929
 930 static float __Covariance(floatVector* function,floatVector* MedianCurve,long i,long around){
 931
 932         long j;
 933         float tmp;
 934
 935         tmp = 0.0;
 936         for(j=(i-around);j<=(i+around);j++){
 937                 tmp += SQR(MedianCurve->data[j] - function->data[j]);
 938         }
 939
 940         return(tmp);
 941 }
 942
 943
 944 void __CovarianceEstimation(floatVector* function,floatVector* MedianCurve,floatVector* Covariance,long around){
 945
 946         long i;
 947         long aroundtmp;
 948
 949         aroundtmp = around;
 950
 951         for(i=0;i<function->size;i++){
 952                 aroundtmp = around;
 953                 if(i!=0 && i!=function->size-1){
 954                         if(i<around){
 955                                 aroundtmp = i;
 956                         }
 957                         else if(i>function->size-1-around){
 958                                 aroundtmp = function->size-1-i;
 959                         }
 960
 961                         Covariance->data[i]  = __Covariance(function,MedianCurve,i,aroundtmp);
 962
 963                 }
 964
 965                 else{
 966                         Covariance->data[i] = 0.0;
 967                 }
 968         }
 969
 970 }
 971
 972
 973 void
 974 lctfDeterminationbyCovarianceMethods(ctfInfo* res, mrcImage* mrc, ctfInfo* ini,long* MaxThonRing, long mode)
 975 {
 976         floatVector* scatter;
 977         floatVector* spacing;
 978         floatVector* ctf;
 979         floatVector* MedianCurve;
 980         floatVector* Covariance;
 981         long         i;
 982         long         Differentialaround;
 983         long         Smootharound;
 984         long         CovarianceAround;
 985         long         ThonRing;
 986         float        E;
 987         float        minE = 0.0;
 988         ctfInfo      tmp;
 989
 990         Differentialaround = 1;
 991         Smootharound       = 2;
 992         CovarianceAround   = 2;
 993
 994         scatter = lmrcFSInfoScatteringAngularDistributionAverage(mrc);
 995         spacing = lmrcFSInfoSpacing(mrc);
 996         DEBUGPRINT2("%f %f\n", scatter->data[1], spacing->data[1]);
 997         ctf     = floatVectorInit(NULL, spacing->size);
 998         tmp     = *ini;
 999         *res    = *ini;
1000
1001         MedianCurve = floatVectorInit(NULL,spacing->size);
1002         Covariance  = floatVectorInit(NULL,spacing->size);
1003
1004         fprintf(stdout,"\nCovarianceMethod\n");
1005         if(*MaxThonRing == -1){
1006                 fprintf(stdout,"all ThonRing search\n");
1007         }
1008         else{
1009                 fprintf(stdout,"only ThonRing = %3ld\n",*MaxThonRing);
1010         }
1011
1012         for(tmp.defocus = 0.0; tmp.defocus<2*ini->defocus; tmp.defocus+=ini->defocus/1000.0) {
1013                 for(i=0;i<spacing->size;i++){
1014                         ctf->data[i] = ctfFunction(&tmp, spacing->data[i],mode);
1015                 }
1016                 __Smoothing(scatter,MedianCurve,Smootharound,1);
1017                 __CovarianceEstimation(scatter,MedianCurve,Covariance,CovarianceAround);
1018                 DEBUGPRINT1("%f\n", ctf->data[1]);
1019                 E= 0.0;
1020                 ThonRing = 0;
1021                 for(i=(0+Differentialaround); i<(spacing->size-Differentialaround); i++) {
1022                         if(tmp.CutLow < spacing->data[i]
1023                         &&spacing->data[i] < tmp.CutHigh ) {
1024                                 if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
1025                                         ThonRing += 1;
1026                                         if(ThonRing <= 3){
1027                                                 E += log10(Covariance->data[i] + 1.0);
1028                                         }
1029                                 }
1030                         }
1031                 }
1032
1033                 if(ThonRing != 0){
1034                         //E = E/(float)ThonRing;
1035                         if((minE > E || minE == 0.0) && (*MaxThonRing == ThonRing || *MaxThonRing == -1) ){
1036                                 fprintf(stdout, "minE: %10.6e  %10.6e %10.6e\n", tmp.defocus,E,minE);
1037                                 minE = E;
1038                                 res->defocus         = tmp.defocus;
1039                         }
1040                 }
1041                 fprintf(stdout, "%10.6e %3ld  %10.6e \n", tmp.defocus,ThonRing, E);
1042         }
1043         fprintf(stdout,"defocus= %10.6f\n",res->defocus);
1044
1045         /*for smooth test*/
1046         for(i=0;i<spacing->size;i++){
1047                 ctf->data[i]  = ctfFunction(res, spacing->data[i],mode);
1048         }
1049         fprintf(stdout,"\n\ndefocus:%f\n",res->defocus);
1050         fprintf(stdout,"    spacing:    SQR(ctf):     scatter:      Median:  Estimation:\n");
1051         for(i=0;i<spacing->size;i++){
1052
1053                 if(tmp.CutLow < spacing->data[i]
1054                 &&spacing->data[i] < tmp.CutHigh ) {
1055                         fprintf(stdout," %10.6e %10.6e %10.6e %10.6e %10.6e ",spacing->data[i],SQR(ctf->data[i]),scatter->data[i],MedianCurve->data[i],Covariance->data[i]);
1056                         if(__SquareLimitValue(ctf,i,Differentialaround) == -1){
1057                                 fprintf(stdout,"ring\n");
1058                         }
1059                         else{
1060                                 fprintf(stdout,"\n");
1061                         }
1062                 }
1063         }
1064 }
1065
1066
1067 void
1068 lctfDeterminationbyMixMethods(ctfInfo* res, mrcImage* mrc, ctfInfo* ini,long* ThonRing, long mode){
1069
1070         fprintf(stdout,"\nDifferential & Covariance Method");
1071         lctfDeterminationbyDifferentialMethods(res, mrc, ini,ThonRing, mode);
1072
1073         lctfDeterminationbyCovarianceMethods(res,mrc,ini,ThonRing,mode);
1074 }