From 048ebbce58111d851dd9cfa7ee8d1e3935e96cff Mon Sep 17 00:00:00 2001
From: Pilar Cossio <pcos@phys01.bc.rzg.mpg.de>
Date: Wed, 11 Nov 2015 11:35:21 +0100
Subject: [PATCH] Polished creat map & others

---
 bioem.cpp       | 50 ++++++++++++++++++++++++++++++++++---------------
 include/param.h |  3 ++-
 model.cpp       |  4 ++--
 param.cpp       | 41 +++++++++++++++++++++++++++++++++++-----
 4 files changed, 75 insertions(+), 23 deletions(-)

diff --git a/bioem.cpp b/bioem.cpp
index 9b3d9dc..1b788a5 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -297,10 +297,12 @@ int bioem::configure(int ac, char* av[])
     {
       if (DebugOutput >= 2 && mpi_rank == 0) timer.ResetStart();
       MPI_Bcast(&param, sizeof(param), MPI_BYTE, 0, MPI_COMM_WORLD);
+      //We have to reinitialize all pointers !!!!!!!!!!!!!!
+      if (mpi_rank != 0) param.angprior = NULL;
 
 //  cout << "HERE MPI " << param.nTotGridAngles << " " << param.nTotGridAngles  * sizeof(myfloat3_t*) << "\n";
-   if (mpi_rank != 0)param.angles =  (myfloat3_t*) mallocchk(3 * param.nTotGridAngles  * sizeof (myfloat3_t*));
-     MPI_Bcast(param.angles, 3 * param.nTotGridAngles  * sizeof (myfloat3_t*),MPI_BYTE, 0, MPI_COMM_WORLD);
+   if (mpi_rank != 0)param.angles =  (myfloat3_t*) mallocchk(param.nTotGridAngles  * sizeof (myfloat3_t));
+     MPI_Bcast(param.angles, param.nTotGridAngles  * sizeof (myfloat3_t),MPI_BYTE, 0, MPI_COMM_WORLD);
 
      //Allocating memory for prior	
 //   if (mpi_rank != 0){//param.angprior = (myfloat_t*)  mallocchk(3* param.nTotGridAngles * sizeof (myfloat_t*));
@@ -698,6 +700,8 @@ int bioem::run()
           if(!param.doquater){ angProbfile <<" RefMap:  MapNumber ; alpha[rad] - beta[rad] - gamma[rad] - logP - cal log Probability + Constant: Numerical Const.+ log (volume) + prior ang\n" ;}
 	  else { angProbfile <<" RefMap:  MapNumber ; q1 - q2 -q3 - logP- cal log Probability + Constant: Numerical Const. + log (volume) + prior ang\n" ;};
 	  angProbfile <<"************************* HEADER:: NOTATION *******************************************\n";
+//          angProbfile <<"Model Used: " << modelfile.c_str() << "\n";
+//          angProbfile <<"Input Used: " << infile.c_str() << "\n";
 	}
       // Output for Cross Correlation File
       ofstream ccProbfile;
@@ -1072,14 +1076,27 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
       for(int n = 0; n < Model.nPointsModel; n++)
 	{
 	  if(Model.points[n].radius < param.pixelSize){
-	    cout << "Error: Radius less than Pixel size: use keyword NO_PROJECT_RADIUS in inputfile\n";
-	    exit(1);
-	  }
+	 //   cout << "Radius less than Pixel size: use keyword NO_PROJECT_RADIUS in inputfile\n";
+          i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
+          j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
+
+          if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
+            {
+              if (DebugOutput >= 0) cout << "WARNING:::: Model Point out of Projection map: " << i << ", " << j << "\n";
+              //              continue;
+              exit(1);
+            }
+
+          localproj[i * param.param_device.NumberPixels + j] += Model.points[n].density;//Model.NormDen;
+
+	 // exit(1);
+	  }else{
+
 	  //Getting Centers of Sphere
-	  i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
-	  j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f);
+	  i = floor(RotatedPointsModel[n].pos[0] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f) -param.shiftX;
+	  j = floor(RotatedPointsModel[n].pos[1] / param.pixelSize + (myfloat_t) param.param_device.NumberPixels / 2.0f + 0.5f) -param.shiftY;
 	  //Getting the radius
-	  irad=int( Model.points[n].radius / param.pixelSize );
+	  irad=int( Model.points[n].radius / param.pixelSize ) + 1;
 	  rad2= Model.points[n].radius * Model.points[n].radius;
 
 	  if (i < 0 || j < 0 || i >= param.param_device.NumberPixels || j >= param.param_device.NumberPixels)
@@ -1087,7 +1104,7 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
               if (DebugOutput >= 0) cout << "WARNING::: Model Point out of Projection map: " << i << ", " << j << "\n";
 	      cout << "Model point " << n << "Rotation: " << iMap <<" "<< RotatedPointsModel[n].pos[0] << " " << RotatedPointsModel[n].pos[1] << " " <<  RotatedPointsModel[n].pos[2] << "\n";
 	      cout << "Original coor " << n <<" " << Model.points[n].point.pos[0] << " " << Model.points[n].point.pos[1] << " " <<Model.points[n].point.pos[2] << "\n";
-	      cout << "WARNING: Angle orient " << n << " " <<  param.angles[iMap].pos[0] << " " <<  param.angles[iMap].pos[1] << " " << param.angles[iMap].pos[2] << "\n";
+	      cout << "WARNING: Angle orient " << n << " " <<  param.angles[iMap].pos[0] << " " <<  param.angles[iMap].pos[1] << " " << param.angles[iMap].pos[2] <<  " out " << i << " " << j << "\n";
 	      cout << "WARNING: MPI rank " << mpi_rank <<"\n";
 	      //              continue;
 	      exit(1);
@@ -1095,20 +1112,23 @@ int bioem::createProjection(int iMap, mycomplex_t* mapFFT)
 
 
 	  //Projecting over the radius
-	  for(int ii= i - irad; ii < i + irad; ii++)
+	  for(int ii= i - irad; ii < i + irad + 1 ; ii++)
 	    {	
-	      for(int jj = j - irad; jj < j + irad; jj++)
+	      for(int jj = j - irad; jj < j + irad + 1 ; jj++)
 		{
-		  dist= ( (myfloat_t) (ii-i)*(ii-i)+(jj-j)*(jj-j) ) *  param.pixelSize ;
+		  dist= ( (myfloat_t) (ii-i)*(ii-i)+(jj-j)*(jj-j) ) *  param.pixelSize *  param.pixelSize ; //at pixel center
 		  if( dist < rad2 )
 		    {
-		      localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
-			//							/ Model.NormDen * 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 	
-			* 3 / (4 * M_PI * rad2 *  Model.points[n].radius); 
+		      localproj[ii * param.param_device.NumberPixels + jj] += param.pixelSize * param.pixelSize * 2 * sqrt( rad2 - dist ) * Model.points[n].density
+			* 3 / (4 * M_PI * Model.points[n].radius * rad2 );
+//			cout << "TT " << ii << " " << jj << " " << localproj[ii * param.param_device.NumberPixels + jj] << " " <<  sqrt( rad2 - dist )  << " " <<  (4 * M_PI * Model.points[n].radius * rad2 )/3. << " " << param.pixelSize *  param.pixelSize * 2 * sqrt( rad2 - dist ) <<  "\n"; 
+                       // This is the number of electrons that this piece of volume contributes (from the full circle)
 		    }
 		}
 	    }
 	}
+
+     }
     }
   else 
     {
diff --git a/include/param.h b/include/param.h
index f0711cb..2f37607 100644
--- a/include/param.h
+++ b/include/param.h
@@ -116,7 +116,8 @@ public:
 	int nTotGridAngles;
 	int nTotCTFs;
 	int nTotCC;
-
+	int shiftX,shiftY;
+	
 	bool printModel;
 	int printModelOrientation;
 	int printModelConvolution;
diff --git a/model.cpp b/model.cpp
index 5ace0cf..9e72375 100644
--- a/model.cpp
+++ b/model.cpp
@@ -222,8 +222,8 @@ int bioem_model::readModel(bioem_param& param, const char* filemodel)
 		r_cm.pos[1] += points[n].point.pos[1]*points[n].density;
 		r_cm.pos[2] += points[n].point.pos[2]*points[n].density;
 	}
-	r_cm.pos[0] = r_cm.pos[0] / NormDen;
-	r_cm.pos[1] = r_cm.pos[1] / NormDen;
+	r_cm.pos[0] = r_cm.pos[0] / NormDen ; 
+	r_cm.pos[1] = r_cm.pos[1] / NormDen ; 
 	r_cm.pos[2] = r_cm.pos[2] / NormDen;
 
 	for(int n = 0; n < nPointsModel; n++)
diff --git a/param.cpp b/param.cpp
index b5b342d..9a06798 100644
--- a/param.cpp
+++ b/param.cpp
@@ -92,6 +92,9 @@ int bioem_param::readParameters(const char* fileinput)
 //f(notuniformangles)inanglef=std::string(fileangles);
   NotUn_angles=0;
   priorMod=1; //Default
+  shiftX=0;
+  shiftY=0;
+
 
   ifstream input(fileinput);
   if (!input.good())
@@ -372,6 +375,18 @@ int bioem_param::readParameters(const char* fileinput)
 	  yespriorAngles=true;
 	  cout << "READING Priors for Orientations in additonal orientation file\n" ;
 	}
+      else if (strcmp(token, "SHIFT_X") == 0)
+        {
+          token = strtok(NULL, " ");
+          shiftX=atoi(token);
+          cout << "Shifting initial model X by "<< shiftX << "\n" ;
+        }
+	else if (strcmp(token, "SHIFT_Y") == 0)
+        {
+          token = strtok(NULL, " ");
+          shiftY=atoi(token);
+          cout << "Shifting initial model Y by "<< shiftY << "\n" ;
+        }
 
 
     }
@@ -456,6 +471,9 @@ int bioem_param::forprintBest(const char* fileinput)
   printrotmod=false;
   readquatlist=false;
   doaaradius=true;
+  shiftX=0;
+  shiftY=0;
+
 
 
   if (!input.good())
@@ -608,6 +626,19 @@ int bioem_param::forprintBest(const char* fileinput)
 	  printrotmod=true;
 	  cout << "PRINTING out rotatted models (best for debugging)\n";
 	}
+     else if (strcmp(token, "SHIFT_X") == 0)
+        {
+          token = strtok(NULL, " ");
+          shiftX=atoi(token);
+          cout << "Shifting initial model X by "<< shiftX << "\n" ;
+        }
+        else if (strcmp(token, "SHIFT_Y") == 0)
+        {
+          token = strtok(NULL, " ");
+          shiftY=atoi(token);
+          cout << "Shifting initial model Y by "<< shiftY << "\n" ;
+        }
+
 
     }
 
@@ -729,7 +760,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
 
       // Euler Angle Array
 
-      angles = (myfloat3_t*) mallocchk( 3 * angleGridPointsAlpha * angleGridPointsBeta * angleGridPointsAlpha * sizeof(myfloat3_t));
+      angles = (myfloat3_t*) mallocchk(  angleGridPointsAlpha * angleGridPointsBeta * angleGridPointsAlpha * sizeof(myfloat3_t));
 
 
       for (int ialpha = 0; ialpha < angleGridPointsAlpha; ialpha ++)
@@ -783,7 +814,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
      // NotUn_angles=NotUn_angles+1;
 
 
-      angles = (myfloat3_t*) mallocchk( NotUn_angles * 3 * sizeof(myfloat3_t));
+      angles = (myfloat3_t*) mallocchk( NotUn_angles * sizeof(myfloat3_t));
 
       if(yespriorAngles){
 	delete[] angprior;
@@ -904,7 +935,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
 
       nTotGridAngles = n;
 
-      angles = (myfloat3_t*) mallocchk( nTotGridAngles * 3 * sizeof(myfloat3_t));
+      angles = (myfloat3_t*) mallocchk( nTotGridAngles * sizeof(myfloat3_t));
  
       voluang= dgridq * dgridq * dgridq * priorMod;
 
@@ -966,7 +997,7 @@ int bioem_param::CalculateGridsParam(const char* fileangles) //TO DO FOR QUATERN
 	}else{
       cout << "Number of quaternions " << ntotquat << "\n";
 	}
-      angles = (myfloat3_t*) mallocchk( ntotquat * 3 * sizeof(myfloat3_t));
+      angles = (myfloat3_t*) mallocchk( ntotquat * sizeof(myfloat3_t));
       //    delete[] angles;
       //    angles = new myfloat3_t[ ntotquat] ;
 
@@ -1236,7 +1267,7 @@ bioem_param::~bioem_param()
   if (refCTF) delete[] refCTF;
   if (CtfParam) delete[] CtfParam;
   if (angles) free(angles);
-//  if(angprior) delete[] angprior;
+  if(angprior) delete[] angprior;
   refCTF= NULL;
   CtfParam = NULL;
   angles = NULL;
-- 
GitLab