Skip to content

GitLab

Commit 048ebbce authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Polished creat map & others

parent 588ee80d
Hide whitespace changes
Inline Side-by-side
Showing with 75 additions and 23 deletions
bioem.cpp
View file @ 048ebbce
......@@ -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
{
......
include/param.h
View file @ 048ebbce
......@@ -116,7 +116,8 @@ public:
int nTotGridAngles;
int nTotCTFs;
int nTotCC;
int shiftX,shiftY;
bool printModel;
int printModelOrientation;
int printModelConvolution;
......
model.cpp
View file @ 048ebbce
......@@ -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++)
......
param.cpp
View file @ 048ebbce
......@@ -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;
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment