"rebuld"

diff --git a/src/Objects/DataManip/mrcImage/src/lmrcFFTCentralSection.h b/src/Objects/DataManip/mrcImage/src/lmrcFFTCentralSection.h
index 10edaa4..7a074f0 100644 (file)
--- a/src/Objects/DataManip/mrcImage/src/lmrcFFTCentralSection.h
+++ b/src/Objects/DataManip/mrcImage/src/lmrcFFTCentralSection.h
@@ -77,6 +77,7 @@ typedef struct lmrcMultiFFTCentralSectionsCompareInfo{
     int PriorSize;
     double Variat;
     double** Sigma;
+    double SigMin;
     lmrcMultiFFTCentralSectionsCompareInfoOut* Out;
        lmrcMultiFFTCentralSectionsCompareInfoOut* Prior;
 }lmrcMultiFFTCentralSectionsCompareInfo;
@@ -106,13 +107,14 @@ extern void lmrcMultiFFTCentralSectionsCompare(lmrcMultiFFTCentralSectionsCompar
 extern void lmrcMultiFFTCentralSectionsCompareNormalization(lmrcMultiFFTCentralSectionsCompareInfoOut* Out, lmrcMultiFFTCentralSectionsCompareInfo* linfo, int mode);
 extern void lmrcMultiFFTCentralSectionsCompareInfoWrite(FILE* fpt, char* filename, lmrcMultiFFTCentralSectionsCompareInfoOut* Out, lmrcMultiFFTCentralSectionsCompareInfo* linfo, float mode1, float mode2);
 extern void lmrcMultiFFTCentralSectionsCompareInfoSort(lmrcMultiFFTCentralSectionsCompareInfoOut Out[], int left, int right);
-//extern void lmrcMultiFFTCentralSectionsCompareInfoSortPre(lmrcMultiFFTCentralSectionsCompareInfo* linfo, int left, int right);
 extern void lmrcMultiFFTCentralSectionsCompareInfoProbSet(lmrcMultiFFTCentralSectionsCompareInfo* linfo ,int mode); 
 extern void lmrcMultiFFTCentralSectionsCompareInfoLimit(lmrcMultiFFTCentralSectionsCompareInfoOut* Out, lmrcMultiFFTCentralSectionsCompareInfo* linfo, float mode1, float mode2);
 extern void lmrcMultiFFTCentralSectionsCompareInfoUpdate(lmrcMultiFFTCentralSectionsCompareInfoOut* Out, lmrcMultiFFTCentralSectionsCompareInfo* linfo);
 extern void lmrcMultiFFTCentralSectionsCompareInfoVariation(lmrcMultiFFTCentralSectionsCompareInfoOut* Out, lmrcMultiFFTCentralSectionsCompareInfo* linfo);
+extern void lmrcMultiFFTCentralSectionsCompareSigma(mrcImage* in, lmrcMultiFFTCentralSectionsCompareInfo* linfo, int mode);
 /* prototype end */
 
+
 #ifdef __cplusplus
 };
 #endif

diff --git a/src/Objects/DataManip/mrcImage/src/lmrcMultiFFTCentralSectionsCompare.c b/src/Objects/DataManip/mrcImage/src/lmrcMultiFFTCentralSectionsCompare.c
index 5c89298..619775d 100644 (file)
--- a/src/Objects/DataManip/mrcImage/src/lmrcMultiFFTCentralSectionsCompare.c
+++ b/src/Objects/DataManip/mrcImage/src/lmrcMultiFFTCentralSectionsCompare.c
@@ -24,54 +24,21 @@ static char __sccs_id[] = "%Z%lmrcMultiFFTCentralSectionsCompare ver%I%; Date:%D
 void
 lmrcMultiFFTCentralSectionsCompare(lmrcMultiFFTCentralSectionsCompareInfoOut* Out, mrcImage* in, mrcImage* volume, lmrcMultiFFTCentralSectionsCompareInfo* linfo, int mode, int nummode){
     mrcImageParaTypeReal x, y;
-    double likelihood;
+    double likelihood, prelikelihood;
     double likelihoodxre, likelihoodyre, likelihoodxyre;
     double likelihoodxim, likelihoodyim, likelihoodxyim;
     double rein, imin;
+    double inOrigin, inlim;
     double revol, imvol;
     int i,j, flag;
     clock_t start,end;
     int sig_x,sig_y;
-    double sigma;
+    double sigma, sigMax=0, sigMin=100;;
     
     DEBUGPRINT("lmrcMultiFFTCentralSectionsCompare start\n");
-
-    //Sigma init
-    linfo->Sigma = malloc(sizeof(double*)*(in->HeaderN.x));
-    if(linfo->Sigma ==NULL){
-    DEBUGPRINT("malloc error\n");
-    }
-    sig_x = 0;
-    for(x=-in->HeaderN.x/2.0; x< in->HeaderN.x/2.0; x++){
-       linfo->Sigma[sig_x] = malloc(sizeof(double)*in->HeaderN.y);
-       if(linfo->Sigma[sig_x]==NULL){
-           DEBUGPRINT("malloc error\n");
-       }
-       sig_x ++;
-    } 
     
-    //Sigma calculate
-    sig_x=0; 
-    for(x = -in->HeaderN.x/2.0 ; x < in->HeaderN.x/2.0; x++){
-    sig_y=0;
-    for(y = -in->HeaderN.y/2.0 ; y < in->HeaderN.y/2.0; y++){
-        sigma=0;
-        for(i=0; i < linfo->PriorSize; i++){        
-            mrcPixelDataGet(in, x, y, 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
-            mrcPixelDataGet(in, x, y, 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
-            mrcPixelDataGet(&Out[linfo->Prior[i].OriginNum].out, x, y, 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
-            mrcPixelDataGet(&Out[linfo->Prior[i].OriginNum].out, x, y, 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
-            sigma = sigma + linfo->Prior[i].Prior*((rein - revol)*(rein - revol) - (imin - imvol)*(imin - imvol));
-        }
-        linfo->Sigma[sig_x][sig_y] = sigma/2;
- 
-     if(sqrt(linfo->Sigma[sig_x][sig_y]*linfo->Sigma[sig_x][sig_y])  < 1.0){
-        linfo->Sigma[sig_x][sig_y] = 1.0;
-        }
-        DEBUGPRINT3("sigma %d %d %f\n",sig_x, sig_y, linfo->Sigma[sig_x][sig_y]);
-        sig_y++;
-    }
-    sig_x++;
+    if(mode != 1){
+        lmrcMultiFFTCentralSectionsCompareSigma(in, linfo, 0);//0:sigma=1 1:sigmaCalcu
     }
 
     for(i=0; i < linfo->OutSize; i++){
@@ -89,6 +56,7 @@ lmrcMultiFFTCentralSectionsCompare(lmrcMultiFFTCentralSectionsCompareInfoOut* Ou
                 for(y=-in->HeaderN.y/2.0; y < in -> HeaderN.y/2.0 ; y++){
                     
                     mrcPixelDataGet(in, x, y, 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                    if(rein > 0){
                     mrcPixelDataGet(in, x, y, 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
                     mrcPixelDataGet(&Out[i].out, x, y, 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
                     mrcPixelDataGet(&Out[i].out, x, y, 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
@@ -99,35 +67,46 @@ lmrcMultiFFTCentralSectionsCompare(lmrcMultiFFTCentralSectionsCompareInfoOut* Ou
                     likelihoodyim = likelihoodyim + imvol*imvol;
                     likelihoodxyre = likelihoodxyre + rein*revol;
                     likelihoodxyim = likelihoodxyim + imin*imvol;
+                    }
                 }
                 }
                 Out[i].Likelihood = likelihoodxyre/sqrt(likelihoodxre*likelihoodyre) + likelihoodxyim/sqrt(likelihoodxim*likelihoodyim);
             }else{
-                likelihood = 0;
+                likelihood = 1;
+                prelikelihood = 0;
                 sigma = 1;
                 sig_x =0;
                 for(x= -in->HeaderN.x/2.0; x < in -> HeaderN.x/2.0 ; x++){
                     sig_y=0;
                 for(y= -in->HeaderN.y/2.0; y < in -> HeaderN.y/2.0 ; y++){
-                    mrcPixelDataGet(in, x, y, 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
-                    mrcPixelDataGet(in, x, y, 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
-                    mrcPixelDataGet(&Out[i].out, x, y, 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
-                    mrcPixelDataGet(&Out[i].out, x, y, 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
-                
-                    likelihood = likelihood + (rein - revol)*(rein - revol) + (imin - imvol)*(imin - imvol)/((-2)*linfo->Sigma[sig_x][sig_y]);
-                    //sigma = sigma*(linfo->Sigma[sig_x][sig_y]);
-                    sigma = sigma*(linfo->Sigma[sig_x][sig_y])*2*M_PI;
-   //               DEBUGPRINT7(" Likelihood  sigma %d %f %f %d %d    %f %f\n", i, x, y, sig_x, sig_y, likelihood, sigma);
+                   mrcPixelDataGet(in, x, y, 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+ //                   mrcPixelDataGet(in, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                    if(rein > 0 && linfo->Sigma[sig_x][sig_y] > 0.0){
+                   mrcPixelDataGet(in, x, y, 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                   mrcPixelDataGet(&Out[i].out, x, y, 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                   mrcPixelDataGet(&Out[i].out, x, y, 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+ //                   mrcPixelDataGet(in, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+   //                 mrcPixelDataGet(&Out[i].out, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+     //               mrcPixelDataGet(&Out[i].out, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                    
+                    
+                   // sigma = ((rein - revol)*(rein - revol) + (imin - imvol)*(imin - imvol))/((-2)*linfo->Sigma[sig_x][sig_y]);
+                   //// prelikelihood = (linfo->SigMin*2*M_PI)*(1/((linfo->Sigma[sig_x][sig_y])*2*M_PI))*(exp(((rein - revol)*(rein - revol) + (imin - imvol)*(imin - imvol))/((-2)*linfo->Sigma[sig_x][sig_y])));
+                    prelikelihood =  prelikelihood + (((rein - revol)*(rein - revol) + (imin - imvol)*(imin - imvol))/((-2)*linfo->Sigma[sig_x][sig_y]));
+  //                 if(prelikelihood > 0.0){
+    //                likelihood = likelihood * prelikelihood;
+      //             }
+        //            sigma= sigma * linfo->Sigma[sig_x][sig_y];
+                 // DEBUGPRINT5(" preLikelihood likelihood, Sigma %d    %e  %e   %e  %e \n", i, prelikelihood, sigma, likelihood, linfo->Sigma[sig_x][sig_y]);
+                    }
                     sig_y ++;
                 }
                 sig_x ++;
                 }
-                DEBUGPRINT2("sigma likelihood %f %f\n", linfo->Sigma[sig_x-1][sig_y-1], likelihood);
- //               if(sigma < 1.0){
-   //                 sigma = 1.0;
-     //           }
-                //Out[i].Likelihood = exp(likelihood)/2*M_PI*sigma;
-                Out[i].Likelihood = exp(likelihood)/sigma;
+ //               DEBUGPRINT3("sigma likelihood %f %f %f\n", linfo->Sigma[sig_x-1][sig_y-1], prelikelihood, likelihood);
+                  Out[i].Likelihood = exp(prelikelihood)/(2*M_PI);
+               //       Out[i].Likelihood = likelihood/(2*M_PI);
+                //   Out[i].Likelihood = likelihood;
             }
             flag = 1;
             }
@@ -135,7 +114,7 @@ lmrcMultiFFTCentralSectionsCompare(lmrcMultiFFTCentralSectionsCompareInfoOut* Ou
         if(flag == 0){
             Out[i].Likelihood = 0;
         }
-    DEBUGPRINT2("Likelihood  %d %f\n", i, Out[i].Likelihood);
+ //  DEBUGPRINT3("Likelihood  %d     %f   %e\n", i, prelikelihood, Out[i].Likelihood);
     }
     DEBUGPRINT("lmrcMultiFFTCentralSectionsCompare end\n");
 }
@@ -184,7 +163,7 @@ lmrcMultiFFTCentralSectionsCompareNormalization(lmrcMultiFFTCentralSectionsCompa
         for(i=0; i<linfo->OutSize; i++){
             sum = sum + Out[i].Post;
         }
-   //     DEBUGPRINT1("lmrcMultiFFTCentralSectionsComapreNormalization sum: %f\n",sum);
+        DEBUGPRINT1("lmrcMultiFFTCentralSectionsComapreNormalization sum: %f\n",sum);
         for(i=0; i<linfo->OutSize; i++){
             if(sum != 0){
             Out[i].Prob = Out[i].Post / sum;
@@ -299,3 +278,88 @@ lmrcMultiFFTCentralSectionsCompareInfoVariation(lmrcMultiFFTCentralSectionsCompa
     linfo->Variat = var;
     printf("Variation: %f\n",linfo->Variat);
 }
+
+void 
+lmrcMultiFFTCentralSectionsCompareSigma(mrcImage* in, lmrcMultiFFTCentralSectionsCompareInfo* linfo, int mode){
+
+    mrcImageParaTypeReal x, y;
+    double rein, imin;
+    double revol, imvol;
+    int i,j;
+    int sig_x, sig_y;
+    double sigma, sigMax=0, sigMin=100, delta;
+    double inOrigin, inlim;
+
+    linfo->Sigma = malloc(sizeof(double*)*(in->HeaderN.x));
+    if(linfo->Sigma ==NULL){
+        DEBUGPRINT("malloc error\n");
+    }    
+    sig_x = 0;
+    for(x=-in->HeaderN.x/2.0; x< in->HeaderN.x/2.0; x++){
+        linfo->Sigma[sig_x] = malloc(sizeof(double)*in->HeaderN.y);
+        if(linfo->Sigma[sig_x]==NULL){
+            DEBUGPRINT("malloc error\n");
+        } 
+        sig_x ++;
+    }    
+    if(mode == 1){
+        mrcPixelDataGet(in, 0, 0, 0, &inOrigin, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+        inlim = inOrigin / 1000;
+        
+        sig_x=0;
+        for(x = -in->HeaderN.x/2.0 ; x < in->HeaderN.x/2.0; x++){
+        sig_y=0;
+        for(y = -in->HeaderN.y/2.0 ; y < in->HeaderN.y/2.0; y++){
+            sigma=0;
+            mrcPixelDataGet(in, x, y, 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+ //           mrcPixelDataGet(in, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+            if(rein > 0.0){
+                for(i=0; i < linfo->PriorSize; i++){
+                   mrcPixelDataGet(in , x, y, 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                   mrcPixelDataGet(&linfo->Out[i].out, x, y, 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                   mrcPixelDataGet(&linfo->Out[i].out, x, y, 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                //    mrcPixelDataGet(in, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &rein, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                //    mrcPixelDataGet(in, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &imin, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                //    mrcPixelDataGet(&linfo->Out[linfo->Prior[i].OriginNum].out, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &revol, mrcPixelRePart, linfo->llinfo.llinfo.InterpMode);
+                    mrcPixelDataGet(&linfo->Out[linfo->Prior[i].OriginNum].out, (x/(in->HeaderN.x/2.0))*(x/(in->HeaderN.x/2.0)), (y/(in->HeaderN.y/2.0))*(y/(in->HeaderN.y/2.0)), 0, &imvol, mrcPixelImPart, linfo->llinfo.llinfo.InterpMode);
+                    sigma = sigma + linfo->Prior[i].Prior*((rein - revol)*(rein - revol) + (imin - imvol)*(imin - imvol));
+                }
+                linfo->Sigma[sig_x][sig_y] = sigma/2; 
+            
+                if(sigMax < linfo->Sigma[sig_x][sig_y]){
+                    sigMax = linfo->Sigma[sig_x][sig_y];
+                }else if(sigMin > linfo->Sigma[sig_x][sig_y] && linfo->Sigma[sig_x][sig_y] > 0.0 ){
+                //}else if(sigMin > linfo->Sigma[sig_x][sig_y] && linfo->Sigma[sig_x][sig_y] > 0.000001 ){
+                    sigMin = linfo->Sigma[sig_x][sig_y];
+                }//else if(linfo->Sigma[sig_x][sig_y] < 0.000001){
+                 //   linfo->Sigma[sig_x][sig_y] = 0.0;
+               // }
+            
+            }
+            else{
+                linfo->Sigma[sig_x][sig_y] = 0.0;
+            }
+ //       DEBUGPRINT3("sigma %d %d %f\n",sig_x, sig_y, linfo->Sigma[sig_x][sig_y]);
+            sig_y++;
+        }
+        sig_x++;
+        } 
+        DEBUGPRINT2("sigMax sigMin %e %e\n",sigMax, sigMin);
+      linfo->SigMin = sigMin;
+       // delta = sqrt((sigMax-sigMin)*(sigMax-sigMin));
+ /*       for(sig_x =0; sig_x < (int)in->HeaderN.x; sig_x++){
+        for(sig_y =0; sig_y < (int)in->HeaderN.y; sig_y++){
+               linfo->Sigma[sig_x][sig_y] = linfo->Sigma[sig_x][sig_y]/(sigMin*2*M_PI);
+           //    linfo->Sigma[sig_x][sig_y] = linfo->Sigma[sig_x][sig_y] /(sigMin*M_PI);
+            
+        }
+        }*/   
+    }else{
+        for(sig_x =0; sig_x < (int)in->HeaderN.x; sig_x++){
+        for(sig_y =0; sig_y < (int)in->HeaderN.y; sig_y++){
+            linfo->Sigma[sig_x][sig_y] = 1.0;
+        }
+        }
+    }
+
+}

diff --git a/src/Tools/mrcImage/mrcMultiFFTCentralSectionsCompare/src/mrcMultiFFTCentralSectionsCompare.c b/src/Tools/mrcImage/mrcMultiFFTCentralSectionsCompare/src/mrcMultiFFTCentralSectionsCompare.c
index 35f5c84..e3f7bf5 100755 (executable)
--- a/src/Tools/mrcImage/mrcMultiFFTCentralSectionsCompare/src/mrcMultiFFTCentralSectionsCompare.c
+++ b/src/Tools/mrcImage/mrcMultiFFTCentralSectionsCompare/src/mrcMultiFFTCentralSectionsCompare.c
@@ -158,6 +158,10 @@ main(int argc, char* argv[])
 
     DEBUGPRINT1("Prior: %f\n",linfo.Out[0].Prob);
 
+    for(i=0; i< linfo.PriorSize; i++){
+ if(linfo.Out[i].Prob >0.0){
+        DEBUGPRINT1(" %15.6e \n", linfo.Out[i].Prob);
+    }}
     free(linfo.Out);
     free(linfo.Prior);