diff --git a/bioem.cpp b/bioem.cpp index ee6c264e307a27c06b43949cfcb09363844c4e2a..79625bba8260fa1ec477ab67ce083ee1232ffaa9 100755 --- a/bioem.cpp +++ b/bioem.cpp @@ -57,8 +57,8 @@ int bioem::configure(int ac, char* av[]) param.writeAngles = false; param.dumpMap = false; param.loadMap = false; - RefMap.readMRC = false; - RefMap.readMultMRC = false; + RefMap.readMRC = false; + RefMap.readMultMRC = false; // ************************************************************************************* cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n"; @@ -68,17 +68,17 @@ int bioem::configure(int ac, char* av[]) try { po::options_description desc("Command line inputs"); - desc.add_options() - ("Inputfile", po::value(), "(Mandatory) Name of input parameter file") - ("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file") - ("Particlesfile", po::value< std::string>(), "(Mandatory) 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") - ; + desc.add_options() + ("Inputfile", po::value(), "(Mandatory) Name of input parameter file") + ("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file") + ("Particlesfile", po::value< std::string>(), "(Mandatory) 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") + ; po::positional_options_description p; @@ -86,8 +86,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("ReadMRC", -1); + p.add("ReadMultipleMRC", -1); p.add("DumpMaps", -1); p.add("LoadMapDump", -1); @@ -125,17 +125,17 @@ int bioem::configure(int ac, char* av[]) 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("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")) { @@ -379,8 +379,8 @@ int bioem::run() int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap) { - //*************************************************************************************** - //***** BioEM routine for comparing reference maps to convoluted maps ***** + //*************************************************************************************** + //***** BioEM routine for comparing reference maps to convoluted maps ***** if (FFTAlgo) { //With FFT Algorithm @@ -405,8 +405,8 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC) { - //*************************************************************************************** - //***** Calculating cross correlation in FFTALGOrithm ***** + //*************************************************************************************** + //***** Calculating cross correlation in FFTALGOrithm ***** const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize]; for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++) diff --git a/include/map.h b/include/map.h index 796a5e3e3be93150cb6f5ae9f72fd0d307d6aa7e..05a11ed894f01802734daf375619ed4785e7232c 100755 --- a/include/map.h +++ b/include/map.h @@ -34,16 +34,16 @@ public: int precalculate(bioem_param& param, bioem& bio); 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); + 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; - bool readMRC,readMultMRC; + bool readMRC,readMultMRC; int ntotRefMap; int numPixels; diff --git a/map.cpp b/map.cpp index a22e1cd9d15db40bcdb4fafadf2ca1dc57374eb5..789243583c3f3ee5b1f522759ba2dcd2058df889 100755 --- a/map.cpp +++ b/map.cpp @@ -35,54 +35,54 @@ int bioem_RefMap::readRefMaps(bioem_param& param) cout << "Particle Maps read from Map Dump\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 << "\n+++++++++++++++++++++++++++++++++++++++++++ \n"; - cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ; - } else { - read_MRC(filemap,param); - cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n"; - cout << "Particle Maps read from ONE MRC File: " << filemap << "\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 << "\n+++++++++++++++++++++++++++++++++++++++++++ \n"; + cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ; + } else { + read_MRC(filemap,param); + cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n"; + cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ; + } + } else { int nummap = -1; @@ -177,7 +177,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param) } cout << "Total Number of particles: " << ntotRefMap; - cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n"; + cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n"; sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap); sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap); @@ -265,266 +265,266 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio) set_pointers(); } -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+++++++++++++++++++++++++++++++++++++++++++\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 = %4d (only data type mode 2: 32-bit)\n",mode); - printf("NSYMBT = %4d (# 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]; - cptr[3]=tmp; - tmp = cptr[1]; - cptr[1]=cptr[2]; - cptr[2]=tmp; - } - return 1; -} -int bioem_RefMap::read_float_empty (FILE *fin) { - float currfloat; - - if (fread(&currfloat,4,1,fin)!=1) return 0; - return 1; -} - -int bioem_RefMap::read_char_float (float *currfloat, FILE *fin) { - char currchar; - - if (fread(&currchar,1,1,fin)!=1) return 0; - *currfloat=(float)currchar; - return 1; -} - -int bioem_RefMap::test_mrc (const char *vol_file, int swap) { - FILE *fin; - int nc, nr, ns, mx, my, mz; - int mode, ncstart, nrstart, nsstart; - float xlen, ylen, zlen; - int i, header_ok = 1, n_range_viols = 0; - int mapc, mapr, maps; - float alpha, beta, gamma; - float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin; - - fin = fopen(vol_file, "rb"); - if( fin == NULL ) { - cout << "ERROR opening MRC: " << vol_file; - exit(1); - } - - //* read header info - 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(&alpha,fin,swap); - header_ok *= read_float(&beta,fin,swap); - header_ok *= read_float(&gamma,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); - for (i=23; i<50; ++i) header_ok *= read_float(&dummy,fin,swap); - header_ok *= read_float(&xorigin,fin,swap); - header_ok *= read_float(&yorigin,fin,swap); - header_ok *= read_float(&zorigin,fin,swap); - fclose (fin); - if (header_ok == 0) { - cout << "ERROR reading MRC header: " << vol_file; - exit(1); - } - - n_range_viols += (nc>5000); n_range_viols += (nc<0); - n_range_viols += (nr>5000); n_range_viols += (nr<0); - n_range_viols += (ns>5000); n_range_viols += (ns<0); - n_range_viols += (ncstart>5000); n_range_viols += (ncstart<-5000); - n_range_viols += (nrstart>5000); n_range_viols += (nrstart<-5000); - n_range_viols += (nsstart>5000); n_range_viols += (nsstart<-5000); - n_range_viols += (mx>5000); n_range_viols += (mx<0); - n_range_viols += (my>5000); n_range_viols += (my<0); - n_range_viols += (mz>5000); n_range_viols += (mz<0); - n_range_viols += (alpha>360.0f); n_range_viols += (alpha<-360.0f); - n_range_viols += (beta>360.0f); n_range_viols += (beta<-360.0f); - n_range_viols += (gamma>360.0f); n_range_viols += (gamma<-360.0f); - - return n_range_viols; -} +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+++++++++++++++++++++++++++++++++++++++++++\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 = %4d (only data type mode 2: 32-bit)\n",mode); + printf("NSYMBT = %4d (# 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]; + cptr[3]=tmp; + tmp = cptr[1]; + cptr[1]=cptr[2]; + cptr[2]=tmp; + } + return 1; +} +int bioem_RefMap::read_float_empty (FILE *fin) { + float currfloat; + + if (fread(&currfloat,4,1,fin)!=1) return 0; + return 1; +} + +int bioem_RefMap::read_char_float (float *currfloat, FILE *fin) { + char currchar; + + if (fread(&currchar,1,1,fin)!=1) return 0; + *currfloat=(float)currchar; + return 1; +} + +int bioem_RefMap::test_mrc (const char *vol_file, int swap) { + FILE *fin; + int nc, nr, ns, mx, my, mz; + int mode, ncstart, nrstart, nsstart; + float xlen, ylen, zlen; + int i, header_ok = 1, n_range_viols = 0; + int mapc, mapr, maps; + float alpha, beta, gamma; + float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin; + + fin = fopen(vol_file, "rb"); + if( fin == NULL ) { + cout << "ERROR opening MRC: " << vol_file; + exit(1); + } + + //* read header info + 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(&alpha,fin,swap); + header_ok *= read_float(&beta,fin,swap); + header_ok *= read_float(&gamma,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); + for (i=23; i<50; ++i) header_ok *= read_float(&dummy,fin,swap); + header_ok *= read_float(&xorigin,fin,swap); + header_ok *= read_float(&yorigin,fin,swap); + header_ok *= read_float(&zorigin,fin,swap); + fclose (fin); + if (header_ok == 0) { + cout << "ERROR reading MRC header: " << vol_file; + exit(1); + } + + n_range_viols += (nc>5000); n_range_viols += (nc<0); + n_range_viols += (nr>5000); n_range_viols += (nr<0); + n_range_viols += (ns>5000); n_range_viols += (ns<0); + n_range_viols += (ncstart>5000); n_range_viols += (ncstart<-5000); + n_range_viols += (nrstart>5000); n_range_viols += (nrstart<-5000); + n_range_viols += (nsstart>5000); n_range_viols += (nsstart<-5000); + n_range_viols += (mx>5000); n_range_viols += (mx<0); + n_range_viols += (my>5000); n_range_viols += (my<0); + n_range_viols += (mz>5000); n_range_viols += (mz<0); + n_range_viols += (alpha>360.0f); n_range_viols += (alpha<-360.0f); + n_range_viols += (beta>360.0f); n_range_viols += (beta<-360.0f); + n_range_viols += (gamma>360.0f); n_range_viols += (gamma<-360.0f); + + return n_range_viols; +} diff --git a/param.cpp b/param.cpp index 58b544ad5fbc039fd176f40ba14b23a38997836b..6e514d2704c4ece391ebde21e6128a8db9df0517 100755 --- a/param.cpp +++ b/param.cpp @@ -42,22 +42,22 @@ int bioem_param::readParameters() // ***************************** Reading Input Parameters ****************************** // ************************************************************************************** - // Control for Parameters - bool yesPixSi = false; - bool yesNumPix = false; - bool yesGPal = false; - bool yesGPbe = false; - bool yesGPEnv = false; - bool yesGPamp = false; - bool yesGPpha = false; - bool yesSTEnv = false; - bool yesSTamp = false; - bool yesSTpha = false; - bool yesGSPamp = false ; - bool yesGSPEnv = false ; - bool yesGSPpha = false ; - bool yesMDC = false ; - bool yesGCen = false ; + // Control for Parameters + bool yesPixSi = false; + bool yesNumPix = false; + bool yesGPal = false; + bool yesGPbe = false; + bool yesGPEnv = false; + bool yesGPamp = false; + bool yesGPpha = false; + bool yesSTEnv = false; + bool yesSTamp = false; + bool yesSTpha = false; + bool yesGSPamp = false ; + bool yesGSPEnv = false ; + bool yesGSPpha = false ; + bool yesMDC = false ; + bool yesGCen = false ; ifstream input(fileinput); if (!input.good()) @@ -87,116 +87,116 @@ int bioem_param::readParameters() { token = strtok(NULL, " "); pixelSize = atof(token); - if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);} - cout << "Pixel Sixe " << pixelSize << "\n"; - yesPixSi= true; + if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);} + cout << "Pixel Sixe " << pixelSize << "\n"; + yesPixSi= true; } else if (strcmp(token, "NUMBER_PIXELS") == 0) { token = strtok(NULL, " "); param_device.NumberPixels = int(atoi(token)); - if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);} - cout << "Number of Pixels " << param_device.NumberPixels << "\n"; - yesNumPix= true ; + if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);} + cout << "Number of Pixels " << param_device.NumberPixels << "\n"; + yesNumPix= true ; } else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0) { token = strtok(NULL, " "); angleGridPointsAlpha = int(atoi(token)); - if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);} - cout << "Grid points alpha " << angleGridPointsAlpha << "\n"; - yesGPal= true; + if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);} + cout << "Grid points alpha " << angleGridPointsAlpha << "\n"; + yesGPal= true; } else if (strcmp(token, "GRIDPOINTS_BETA") == 0) { token = strtok(NULL, " "); angleGridPointsBeta = int(atoi(token)); - if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);} - cout << "Grid points beta " << angleGridPointsBeta << "\n"; - yesGPbe= true; + if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);} + cout << "Grid points beta " << angleGridPointsBeta << "\n"; + yesGPbe= true; } else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0) { token = strtok(NULL, " "); numberGridPointsEnvelop = int(atoi(token)); - if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);} - cout << "Grid points envelope " << numberGridPointsEnvelop << "\n"; - yesGPEnv = true; + if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);} + cout << "Grid points envelope " << numberGridPointsEnvelop << "\n"; + yesGPEnv = true; } else if (strcmp(token, "START_ENVELOPE") == 0) { token = strtok(NULL, " "); startGridEnvelop = atof(token); - if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START_ENVELOPE "; exit(1);} - cout << "Start Envelope " << startGridEnvelop << "\n"; - yesSTEnv = true ; + if (startGridEnvelop < 0 ) { cout << "*** Error: Negative START_ENVELOPE "; exit(1);} + cout << "Start Envelope " << startGridEnvelop << "\n"; + yesSTEnv = true ; } else if (strcmp(token, "GRIDSPACE_ENVELOPE") == 0) { token = strtok(NULL, " "); gridEnvelop = atof(token); - if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);} - cout << "Grid spacing Envelope " << gridEnvelop << "\n"; - yesGSPEnv = true ; + if (gridEnvelop < 0 ) { cout << "*** Error: Negative GRIDSPACE_ENVELOPE "; exit(1);} + cout << "Grid spacing Envelope " << gridEnvelop << "\n"; + yesGSPEnv = true ; + } + else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0) + { + token = strtok(NULL," "); + numberGridPointsCTF_phase=int(atoi(token)); + cout << "Grid points PSF " << numberGridPointsCTF_phase << "\n"; + yesGPpha = true; + } + else if (strcmp(token,"START_PSF_PHASE")==0) + { + token = strtok(NULL," "); + startGridCTF_phase=atof(token); + cout << "Start PSF " << startGridCTF_phase << "\n"; + yesSTpha = true ; + } + else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0) + { + token = strtok(NULL," "); + gridCTF_phase=atof(token); + cout << "Grid Space PSF " << gridCTF_phase << "\n"; + yesGSPpha = true ; + } + else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0) + { + token = strtok(NULL," "); + numberGridPointsCTF_amp=int(atoi(token)); + cout << "Grid points PSF " << numberGridPointsCTF_amp << "\n"; + yesGPamp = true ; + } + else if (strcmp(token,"START_PSF_AMP")==0) + { + token = strtok(NULL," "); + startGridCTF_amp=atof(token); + cout << "Start PSF " << startGridCTF_amp << "\n"; + yesSTamp = true ; + } + else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0) + { + token = strtok(NULL," "); + gridCTF_amp=atof(token); + cout << "Grid Space PSF " << gridCTF_amp << "\n"; + yesGSPamp = true ; } - else if (strcmp(token,"GRIDPOINTS_PSF_PHASE")==0) - { - token = strtok(NULL," "); - numberGridPointsCTF_phase=int(atoi(token)); - cout << "Grid points PSF " << numberGridPointsCTF_phase << "\n"; - yesGPpha = true; - } - else if (strcmp(token,"START_PSF_PHASE")==0) - { - token = strtok(NULL," "); - startGridCTF_phase=atof(token); - cout << "Start PSF " << startGridCTF_phase << "\n"; - yesSTpha = true ; - } - else if (strcmp(token,"GRIDSPACE_PSF_PHASE")==0) - { - token = strtok(NULL," "); - gridCTF_phase=atof(token); - cout << "Grid Space PSF " << gridCTF_phase << "\n"; - yesGSPpha = true ; - } - else if (strcmp(token,"GRIDPOINTS_PSF_AMP")==0) - { - token = strtok(NULL," "); - numberGridPointsCTF_amp=int(atoi(token)); - cout << "Grid points PSF " << numberGridPointsCTF_amp << "\n"; - yesGPamp = true ; - } - else if (strcmp(token,"START_PSF_AMP")==0) - { - token = strtok(NULL," "); - startGridCTF_amp=atof(token); - cout << "Start PSF " << startGridCTF_amp << "\n"; - yesSTamp = true ; - } - else if (strcmp(token,"GRIDSPACE_PSF_AMP")==0) - { - token = strtok(NULL," "); - gridCTF_amp=atof(token); - cout << "Grid Space PSF " << gridCTF_amp << "\n"; - yesGSPamp = true ; - } else if (strcmp(token, "MAX_D_CENTER") == 0) { token = strtok(NULL, " "); param_device.maxDisplaceCenter = int(atoi(token)); - if (param_device.maxDisplaceCenter < 0 ) { cout << "*** Error: Negative MAX_D_CENTER "; exit(1);} - cout << "Maximum displacement Center " << param_device.maxDisplaceCenter << "\n"; - yesMDC = true; + if (param_device.maxDisplaceCenter < 0 ) { cout << "*** Error: Negative MAX_D_CENTER "; exit(1);} + cout << "Maximum displacement Center " << param_device.maxDisplaceCenter << "\n"; + yesMDC = true; } else if (strcmp(token, "PIXEL_GRID_CENTER") == 0) { token = strtok(NULL, " "); param_device.GridSpaceCenter = int(atoi(token)); - if (param_device.GridSpaceCenter < 0 ) { cout << "*** Error: Negative PIXEL_GRID_CENTER "; exit(1);} - cout << "Grid space displacement center " << param_device.GridSpaceCenter << "\n"; - yesGCen = true; + if (param_device.GridSpaceCenter < 0 ) { cout << "*** Error: Negative PIXEL_GRID_CENTER "; exit(1);} + cout << "Grid space displacement center " << param_device.GridSpaceCenter << "\n"; + yesGCen = true; } else if (strcmp(token, "WRITE_PROB_ANGLES") == 0) { @@ -206,41 +206,41 @@ int bioem_param::readParameters() } input.close(); - // Checks for ALL INPUT - if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);}; - if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);}; - if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);}; - if( not ( yesGPbe ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);}; - if( not ( yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);}; - if( not ( yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);}; - if( not ( yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);}; - if( not ( yesSTEnv ) ) { cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ; exit (1);}; - if( not ( yesSTamp ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ; exit (1);}; - if( not ( yesSTpha ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ; exit (1);}; - if( not ( yesGSPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ; exit (1);}; - if( not ( yesGSPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ; exit (1);}; - if( not ( yesGSPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ; exit (1);}; - if( not ( yesMDC ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);}; - if( not ( yesGCen ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);}; - - //More checks with input parameters - - // Envelope should not have a standard deviation greater than Npix/2 - if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) { - cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length " ; - exit(1); - } - // Envelop param should be positive - if( startGridCTF_amp < 0 || startGridCTF_amp > 1){ - cout << "PSF Amplitud should be between 0 and 1 \n" ; - exit(1); - } - - if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1){ - cout << "PSF Amplitud should be between 0 and 1 \n" ; - exit(1); - } - + // Checks for ALL INPUT + if( not ( yesPixSi ) ){ cout << "**** INPUT MISSING: Please provide PIXEL_SIZE\n" ; exit (1);}; + if( not ( yesNumPix ) ){ cout << "**** INPUT MISSING: Please provide NUMBER_PIXELS \n" ; exit (1);}; + if( not ( yesGPal ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ALPHA \n" ; exit (1);}; + if( not ( yesGPbe ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_BETA \n" ; exit (1);}; + if( not ( yesGPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_ENVELOPE \n" ; exit (1);}; + if( not ( yesGPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_AMP \n" ; exit (1);}; + if( not ( yesGPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDPOINTS_PSF_PHASE \n" ; exit (1);}; + if( not ( yesSTEnv ) ) { cout << "**** INPUT MISSING: Please provide START_ENVELOPE \n" ; exit (1);}; + if( not ( yesSTamp ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_AMP \n" ; exit (1);}; + if( not ( yesSTpha ) ) { cout << "**** INPUT MISSING: Please provide START_PSF_PHASE \n" ; exit (1);}; + if( not ( yesGSPamp ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_AMP \n" ; exit (1);}; + if( not ( yesGSPEnv ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_ENVELOPE \n" ; exit (1);}; + if( not ( yesGSPpha ) ) { cout << "**** INPUT MISSING: Please provide GRIDSPACE_PSF_PHASE \n" ; exit (1);}; + if( not ( yesMDC ) ) { cout << "**** INPUT MISSING: Please provide MAX_D_CENTER \n" ; exit (1);}; + if( not ( yesGCen ) ) { cout << "**** INPUT MISSING: Please provide PIXEL_GRID_CENTER \n" ; exit (1);}; + + //More checks with input parameters + + // Envelope should not have a standard deviation greater than Npix/2 + if(sqrt(1./( (myfloat_t) numberGridPointsDisplaceCenter * gridEnvelop + startGridEnvelop))> float(param_device.NumberPixels)/2.0) { + cout << "MAX Standar deviation of envelope is larger than Allowed KERNEL Length " ; + exit(1); + } + // Envelop param should be positive + if( startGridCTF_amp < 0 || startGridCTF_amp > 1){ + cout << "PSF Amplitud should be between 0 and 1 \n" ; + exit(1); + } + + if( (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp < 0 || (myfloat_t) numberGridPointsCTF_amp * gridCTF_amp + startGridCTF_amp > 1){ + cout << "PSF Amplitud should be between 0 and 1 \n" ; + exit(1); + } + cout << " +++++++++++++++++++++++++++++++++++++++++ \n";