Skip to content

GitLab

Commit 92e7d74a authored by Pilar Cossio's avatar Pilar Cossio
Browse files

Modified to read Multiple MRC particle maps

parent 81fab4b2
Hide whitespace changes
Inline Side-by-side
Showing with 486 additions and 140 deletions
CMakeLists.txt
View file @ 92e7d74a
......@@ -11,7 +11,7 @@ FIND_PACKAGE(PkgConfig REQUIRED)
PKG_CHECK_MODULES(FFTW fftw3 REQUIRED)
INCLUDE_DIRECTORIES(${FFTW_INCLUDE_DIRS})
FIND_PACKAGE(Boost 1.54 REQUIRED)
FIND_PACKAGE(Boost 1.43 REQUIRED)
INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS})
FIND_PACKAGE(OpenMP)
......
bioem.cpp
View file @ 92e7d74a
......@@ -59,6 +59,8 @@ int bioem::configure(int ac, char* av[])
param.writeAngles=false;
RefMap.dumpMap = false;
RefMap.loadMap = false;
RefMap.readMRC = false;
RefMap.readMultMRC = false;
/*************************************************************************************/
cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
......@@ -73,6 +75,8 @@ int bioem::configure(int ac, char* av[])
("Modelfile", po::value< std::string>() , "Name of model file")
("Particlesfile", po::value< std::string>(), "Name of paricles file")
("ReadPDB", "(Optional) If reading model file in PDB format")
("ReadMRC", "(Optional) If reading particle file in MRC format")
("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
("DumpMaps", "(Optional) Dump maps after they were red from maps file")
("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
("help", "(Optional) Produce help message")
......@@ -83,6 +87,8 @@ int bioem::configure(int ac, char* av[])
p.add("Modelfile", -1);
p.add("Particlesfile", -1);
p.add("ReadPDB", -1);
p.add("ReadMRC", -1);
p.add("ReadMultipleMRC", -1);
p.add("DumpMaps", -1);
p.add("LoadMapDump", -1);
......@@ -119,6 +125,17 @@ int bioem::configure(int ac, char* av[])
cout << "Reading model file in PDB format.\n";
Model.readPDB=true;
}
if (vm.count("ReadMRC"))
{
cout << "Reading particle file in MRC format.\n";
RefMap.readMRC=true;
}
if (vm.count("ReadMultipleMRC"))
{
cout << "Reading Multiple MRCs.\n";
RefMap.readMultMRC=true;
}
if (vm.count("DumpMaps"))
{
......
include/map.h
View file @ 92e7d74a
......@@ -22,7 +22,13 @@ public:
}
int readRefMaps(bioem_param& param);
int PreCalculateMapsFFT(bioem_param& param);
int read_int(int *currlong, FILE *fin, int swap);
int read_float(float *currfloat, FILE *fin, int swap);
int read_float_empty (FILE *fin);
int read_char_float (float *currfloat, FILE *fin) ;
int test_mrc (const char *vol_file, int swap);
int read_MRC(const char* filename,bioem_param& param);
mycomplex_t* RefMapsFFT;
const char* filemap;
......@@ -32,7 +38,7 @@ public:
myfloat_t sumsquare_RefMap[BIOEM_MAX_MAPS];
myfloat_t ForLogProbfromRef[BIOEM_MAX_MAPS];
bool dumpMap, loadMap;
bool dumpMap, loadMap, readMRC,readMultMRC;
__host__ __device__ inline myfloat_t get(int map, int x, int y) const {return(Ref[map].points[x][y]);}
__host__ __device__ inline const myfloat_t* getp(int map, int x, int y) const {return(&Ref[map].points[x][y]);}
......
main.cpp
View file @ 92e7d74a
......@@ -27,7 +27,7 @@ int main(int argc, char* argv[])
#pragma omp parallel
{
_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
//_MM_SET_DENORMALS_ZERO_MODE(_MM_DENORMALS_ZERO_ON); //Flush denormals to zero in all OpenMP threads
}
HighResTimer timer;
......
map.cpp
View file @ 92e7d74a
......@@ -13,160 +13,483 @@ using namespace std;
int bioem_RefMap::readRefMaps(bioem_param& param)
{
/**************************************************************************************/
/***********************Reading reference Particle Maps************************/
/**************************************************************************************/
if (loadMap)
/**************************************************************************************/
/***********************Reading reference Particle Maps************************/
/**************************************************************************************/
if (loadMap)
{
FILE* fp = fopen("maps.dump", "rb");
if (fp == NULL)
{
FILE* fp = fopen("maps.dump", "rb");
if (fp == NULL)
{
cout << "Error opening dump file\n";
exit(1);
}
fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
if (ntotRefMap > BIOEM_MAX_MAPS)
{
cout << "BIOEM_MAX_MAPS too small, some maps dropped\n";
ntotRefMap = BIOEM_MAX_MAPS;
}
fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
fclose(fp);
cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
cout << "Error opening dump file\n";
exit(1);
}
else
fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
if (ntotRefMap > BIOEM_MAX_MAPS)
{
int nummap=-1;
int lasti=0;
int lastj=0;
ifstream input(filemap);
if (!input.good())
{
cout << "Particle Maps Failed to open file" << endl ;
exit(1);
}
cout << "BIOEM_MAX_MAPS too small, some maps dropped\n";
ntotRefMap = BIOEM_MAX_MAPS;
}
fread(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
fclose(fp);
char line[512] = {' '};
char tmpLine[512] = {' '};
while (!input.eof())
{
input.getline(line,512);
strncpy(tmpLine,line,strlen(line));
char *token = strtok(tmpLine," ");
if (strcmp(token,"PARTICLE")==0) // to count the number of maps
{
nummap++;
if (nummap % 128 == 0)
{
cout << "..." << nummap << "\n";
}
if (nummap == BIOEM_MAX_MAPS)
{
cout << "BIOEM_MAX_MAPS too small\n";
exit(1);
}
if(lasti!=param.param_device.NumberPixels && lastj!=param.param_device.NumberPixels && nummap>0)
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
exit(1);
}
}
else
{
int i,j;
float z;
char tmpVals[36] = {' '};
strncpy (tmpVals,line,8);
sscanf (tmpVals,"%d",&i);
strncpy (tmpVals,line+8,8);
sscanf (tmpVals,"%d",&j);
strncpy (tmpVals,line+16,12);
sscanf (tmpVals,"%f",&z);
//checking for Map limits
if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
{
Ref[nummap].points[i-1][j-1] = (myfloat_t) z;
lasti=i;
lastj=j;
}
else
{
cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
exit(1);
}
}
}
cout << ".";
ntotRefMap=nummap+1;
cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
if (dumpMap)
{
FILE* fp = fopen("maps.dump", "w+b");
if (fp == NULL)
{
cout << "Error opening dump file\n";
exit(1);
}
fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
fclose(fp);
cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
}
else
{
if(readMRC){
ntotRefMap=0;
if(readMultMRC)
{
ifstream input(filemap);
if (!input.good())
{
cout << "Failed to open file contaning MRC names: " << filemap << "\n";
exit(0);
}
char line[512] = {0};
char mapname[100];
char tmpm[10] = {0};
const char* indifile;
while (!input.eof())
{
input.getline(line,512);
char tmpVals[100] = {' '};
strncpy (tmpVals,line,100);
sscanf (tmpVals,"%99c",mapname);
// Check for last line
strncpy (tmpm,mapname,3);
if(strcmp(tmpm,"XXX")!=0){
indifile = mapname;
//Reading Multiple MRC
read_MRC(indifile,param);
}
for(int i=0;i<3;i++)mapname[i] = 'X';
for(int i=3;i<100;i++)mapname[i] = {0};
}
cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
} else {
read_MRC(filemap,param);
cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
}
}
else {
int nummap=-1;
int lasti=0;
int lastj=0;
ifstream input(filemap);
if (!input.good())
{
cout << "Particle Maps Failed to open file" << endl ;
exit(1);
}
char line[512] = {' '};
char tmpLine[512] = {' '};
while (!input.eof())
{
input.getline(line,512);
strncpy(tmpLine,line,strlen(line));
char *token = strtok(tmpLine," ");
if (strcmp(token,"PARTICLE")==0) // to count the number of maps
{
nummap++;
if (nummap % 128 == 0)
{
cout << "..." << nummap << "\n";
}
if (nummap == BIOEM_MAX_MAPS)
{
cout << "BIOEM_MAX_MAPS too small\n";
exit(1);
}
if(lasti!=param.param_device.NumberPixels && lastj!=param.param_device.NumberPixels && nummap>0)
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
exit(1);
}
}
else
{
int i,j;
float z;
char tmpVals[36] = {' '};
strncpy (tmpVals,line,8);
sscanf (tmpVals,"%d",&i);
strncpy (tmpVals,line+8,8);
sscanf (tmpVals,"%d",&j);
strncpy (tmpVals,line+16,12);
sscanf (tmpVals,"%f",&z);
//checking for Map limits
if(i>0 && i-1 <BIOEM_MAP_SIZE_X&& j>0 && j-1<BIOEM_MAP_SIZE_Y && nummap < BIOEM_MAX_MAPS)
{
Ref[nummap].points[i-1][j-1] = (myfloat_t) z;
lasti=i;
lastj=j;
}
else
{
cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
exit(1);
}
}
}
cout << ".";
ntotRefMap=nummap+1;
cout << "Particle Maps read from Standard File\n" ;
}
cout << "Total Number of particles: " << ntotRefMap ;
cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
if (dumpMap)
{
FILE* fp = fopen("maps.dump", "w+b");
if (fp == NULL)
{
cout << "Error opening dump file\n";
exit(1);
}
fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
fwrite(&Ref[0], sizeof(Ref[0]), ntotRefMap, fp);
fclose(fp);
}
}
return(0);
return(0);
}
int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
/**************************************************************************************/
/********** Routine that pre-calculates Kernels for Convolution **********************/
/************************************************************************************/
/**************************************************************************************/
/********** Routine that pre-calculates Kernels for Convolution **********************/
/************************************************************************************/
myfloat_t* localMap;
myfloat_t* localMap;
localMap= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
localMap= (myfloat_t *) myfftw_malloc(sizeof(myfloat_t) * param.param_device.NumberPixels * param.param_device.NumberPixels);
RefMapsFFT = new mycomplex_t[ntotRefMap * param.RefMapSize];
RefMapsFFT = new mycomplex_t[ntotRefMap * param.RefMapSize];
mycomplex_t* localout;
localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
mycomplex_t* localout;
localout= (mycomplex_t *) myfftw_malloc(sizeof(mycomplex_t) * param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D);
for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
{
//Assigning localMap values to padded Map
for(int i=0; i< param.param_device.NumberPixels; i++)
{
//Assigning localMap values to padded Map
for(int i=0; i< param.param_device.NumberPixels; i++)
{
for(int j=0; j< param.param_device.NumberPixels; j++)
{
localMap[i*param.param_device.NumberPixels+j]=Ref[iRefMap].points[i][j];
}
}
for(int j=0; j< param.param_device.NumberPixels; j++)
{
localMap[i*param.param_device.NumberPixels+j]=Ref[iRefMap].points[i][j];
}
}
//Calling FFT_Forward
myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localMap,localout);
//Calling FFT_Forward
myfftw_execute_dft_r2c(param.fft_plan_r2c_forward,localMap,localout);
// Normalizing and saving the Reference CTFs
mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.RefMapSize];
for(int i=0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
{
RefMap[i][0]=localout[i][0];
RefMap[i][1]=localout[i][1];
}
// Normalizing and saving the Reference CTFs
mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.RefMapSize];
for(int i=0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
{
RefMap[i][0]=localout[i][0];
RefMap[i][1]=localout[i][1];
}
}
myfftw_free(localMap);
myfftw_free(localout);
myfftw_free(localMap);
myfftw_free(localout);
return(0);
}
int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
{
myfloat_t st,st2;
unsigned long count;
FILE *fin;
float currfloat;
int nc, nr, ns, swap, header_ok = 1;
float xlen, ylen, zlen;
int mode, ncstart, nrstart, nsstart, ispg, nsymbt, lskflg;
float a_tmp, b_tmp, g_tmp;
int mx, my, mz,mapc, mapr, maps;
float dmin, dmax, dmean;
int n_range_viol0, n_range_viol1;
fin = fopen(filename, "rb");
if( fin == NULL ) {
cout << "ERROR opening MRC: " << filename;
exit(1);
}
n_range_viol0 = test_mrc(filename,0);
n_range_viol1 = test_mrc(filename,1);
if (n_range_viol0 < n_range_viol1) { /* guess endianism */
swap = 0;
if (n_range_viol0 > 0) {
printf(" Warning: %i header field range violations detected in file %s \n", n_range_viol0,filename);
}
} else {
swap = 1;
if (n_range_viol1 > 0) {
printf("Warning: %i header field range violations detected in file %s \n", n_range_viol1,filename);
}
}
printf("++++++++++++++++++++++++++++++++++++++++++++++\n");
printf("Reading Information from MRC: %s \n", filename);
header_ok *= read_int(&nc,fin,swap);
header_ok *= read_int(&nr,fin,swap);
header_ok *= read_int(&ns,fin,swap);
header_ok *= read_int(&mode,fin,swap);
header_ok *= read_int(&ncstart,fin,swap);
header_ok *= read_int(&nrstart,fin,swap);
header_ok *= read_int(&nsstart,fin,swap);
header_ok *= read_int(&mx,fin,swap);
header_ok *= read_int(&my,fin,swap);
header_ok *= read_int(&mz,fin,swap);
header_ok *= read_float(&xlen,fin,swap);
header_ok *= read_float(&ylen,fin,swap);
header_ok *= read_float(&zlen,fin,swap);
header_ok *= read_float(&a_tmp,fin,swap);
header_ok *= read_float(&b_tmp,fin,swap);
header_ok *= read_float(&g_tmp,fin,swap);
header_ok *= read_int(&mapc,fin,swap);
header_ok *= read_int(&mapr,fin,swap);
header_ok *= read_int(&maps,fin,swap);
header_ok *= read_float(&dmin,fin,swap);
header_ok *= read_float(&dmax,fin,swap);
header_ok *= read_float(&dmean,fin,swap);
header_ok *= read_int(&ispg,fin,swap);
header_ok *= read_int(&nsymbt,fin,swap);
header_ok *= read_int(&lskflg,fin,swap);
printf("Number Columns = %8d \n",nc);
printf("Number Rows = %8d \n",nr);
printf("Number Sections = %8d \n",ns);
printf("MODE = %8d (only data type mode 2: 32-bit)\n",mode);
printf("NSYMBT = %8d (# bytes symmetry operators)\n",nsymbt);
/* printf(" NCSTART = %8d (index of first column, counting from 0)\n",ncstart);
printf("> NRSTART = %8d (index of first row, counting from 0)\n",nrstart);
printf(" NSSTART = %8d (index of first section, counting from 0)\n",nsstart);
printf(" MX = %8d (# of X intervals in unit cell)\n",mx);
printf(" MY = %8d (# of Y intervals in unit cell)\n",my);
printf(" MZ = %8d (# of Z intervals in unit cell)\n",mz);
printf(" X length = %8.3f (unit cell dimension)\n",xlen);
printf(" Y length = %8.3f (unit cell dimension)\n",ylen);
printf(" Z length = %8.3f (unit cell dimension)\n",zlen);
printf(" Alpha = %8.3f (unit cell angle)\n",a_tmp);
printf(" Beta = %8.3f (unit cell angle)\n",b_tmp);
printf(" Gamma = %8.3f (unit cell angle)\n",g_tmp);
printf(" MAPC = %8d (columns axis: 1=X,2=Y,3=Z)\n",mapc);
printf(" MAPR = %8d (rows axis: 1=X,2=Y,3=Z)\n",mapr);
printf(" MAPS = %8d (sections axis: 1=X,2=Y,3=Z)\n",maps);
printf(" DMIN = %8.3f (minimum density value - ignored)\n",dmin);
printf(" DMAX = %8.3f (maximum density value - ignored)\n",dmax);
printf(" DMEAN = %8.3f (mean density value - ignored)\n",dmean);
printf(" ISPG = %8d (space group number - ignored)\n",ispg);
printf(" NSYMBT = %8d (# bytes storing symmetry operators)\n",nsymbt);
printf(" LSKFLG = %8d (skew matrix flag: 0:none, 1:follows)\n",lskflg);*/
if (header_ok == 0) {
cout << "ERROR reading MRC header: " << filename;
exit(1);
}
if(nr!=param.param_device.NumberPixels || nc!=param.param_device.NumberPixels )
{
cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << nc << ", j " << nr << ")" << "\n";
exit(1);
}
if(mode!=2)
{
cout << "ERROR with MRC mode " << mode << "\n";
cout << "Currently mode 2 is the only one allowed" << "\n";
exit(1);
} else
{
rewind (fin);
for (count=0; count<256; ++count) if (read_float_empty(fin)==0) {
cout << "ERROR Converting Data: " << filename;
exit(1);
}
for (count=0; count<(unsigned long)nsymbt; ++count) if (read_char_float(&currfloat,fin)==0) {
cout << "ERROR Converting Data: " << filename;
exit(1);
}
for ( int nmap = 0 ; nmap < ns ; nmap ++ )
{
st=0.0;
st2=0.0;
for ( int j = 0 ; j < nr ; j ++ )
for ( int i = 0 ; i < nc ; i ++ )
{
if (read_float(&currfloat,fin,swap)==0) {
cout << "ERROR Converting Data: " << filename;
exit(1);
} else {
Ref[nmap+ntotRefMap].points[i][j] = currfloat;
st += currfloat;
st2 += currfloat*currfloat;
}
}
//Normaling maps to zero mean and unit standard deviation
st /= float(nr*nc);
st2 = sqrt(st2/float(nr*nc)-st*st);
for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
Ref[nmap+ntotRefMap].points[i][j] = Ref[nmap+ntotRefMap].points[i][j]/st2 - st/st2;
// if(nmap+ntotRefMap==300)cout << i << " " << j << " " << nmap+ntotRefMap << " " << Ref[nmap+ntotRefMap].points[i][j] << "\n";
}
}
ntotRefMap += ns ;
// cout << ntotRefMap << "\n";
}
fclose (fin);
return(0);
}
int bioem_RefMap::read_float(float *currfloat, FILE *fin, int swap) {
unsigned char *cptr, tmp;
if (fread(currfloat,4,1,fin)!=1) return 0;
if (swap == 1) {
cptr = (unsigned char *)currfloat;
tmp = cptr[0];
cptr[0]=cptr[3];
cptr[3]=tmp;
tmp = cptr[1];
cptr[1]=cptr[2];
cptr[2]=tmp;
}
return 1;
}
int bioem_RefMap::read_int(int *currlong, FILE *fin, int swap) {
unsigned char *cptr, tmp;
if (fread(currlong,4,1,fin)!=1) return 0;
if (swap == 1) {
cptr = (unsigned char *)currlong;
tmp = cptr[0];
cptr[0]=cptr[3];