diff --git a/map.cpp b/map.cpp index 0262d8e70752236e0d0ca055ac79954137fdd18a..a9f75dbd9cd2ad2119113dc4f2c67dd3ade475d7 100644 --- a/map.cpp +++ b/map.cpp @@ -86,7 +86,15 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) if(strcmp(tmpm,"XXX")!=0) { indifile=strline.c_str(); - cout << "FILE Name: " << indifile << "\n"; + + size_t foundpos= strline.find(".mrc"); + size_t endpos = strline.find_last_not_of(" \t"); + + if(foundpos > endpos){ + cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n"; + cout << "Warining:::: Are you sure you want to read an MRC? \n"; + } + //Reading Multiple MRC read_MRC(indifile,param); } @@ -99,6 +107,17 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) } else { + + string strfilename(filemap); + + size_t foundpos= strfilename.find(".mrc"); + size_t endpos = strfilename.find_last_not_of(" \t"); + + if(foundpos > endpos){ + cout << "Warining:::: .mrc extension NOT dectected in file name::" << filemap <<" \n"; + cout << "Warining:::: Are you sure you want to read an MRC? \n"; + } + read_MRC(filemap,param); cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n"; cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ; @@ -120,6 +139,8 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) char tmpLine[512] = {0}; bool first=true; + int countpix; + while (!input.eof()) { input.getline(line, 511); @@ -135,9 +156,11 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) first=false; } + if (strcmp(token, "PARTICLE") == 0) // to count the number of maps { nummap++; + countpix=0; if (allocsize == 0) { allocsize = 64; @@ -162,7 +185,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) { int i, j; float z; - + char tmpVals[36] = {0}; strncpy (tmpVals, line, 8); @@ -171,14 +194,16 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) strncpy (tmpVals, line + 8, 8); sscanf (tmpVals, "%d", &j); - strncpy (tmpVals, line + 16, 12); + strncpy (tmpVals, line + 16, 16); sscanf (tmpVals, "%f", &z); //checking for Map limits - if(i > 0 && i - 1 < param.param_device.NumberPixels && j > 0 && j - 1 < param.param_device.NumberPixels) + if(i > -1 && i < param.param_device.NumberPixels && j > -1 && j < param.param_device.NumberPixels) { - maps[nummap * refMapSize + (i - 1) * numPixels + (j - 1)] = (myfloat_t) z; + countpix++; + maps[nummap * refMapSize + i * numPixels + j] = (myfloat_t) z; lasti = i; lastj = j; + // cout << countpix << " " << param.param_device.NumberPixels*param.param_device.NumberPixels << "\n"; } else { @@ -187,7 +212,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) } } } - if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels ) + if(lasti != param.param_device.NumberPixels-1 || lastj != param.param_device.NumberPixels-1 || countpix != param.param_device.NumberPixels*param.param_device.NumberPixels +1 ) { cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n"; exit(1); @@ -197,20 +222,23 @@ int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap) maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap); cout << "Particle Maps read from Standard File \n"; - if (param.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(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp); - fclose(fp); - } } +//Dumping Maps + if (param.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(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp); + fclose(fp); + } + + if (getenv("BIOEM_DEBUG_NMAPS")) { ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS")); @@ -299,7 +327,8 @@ void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio) nAngles = angles; nCC = cc; ptr = bio.malloc_device_host(get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)); - cout << "Allocation #Maps " << maps << " #Angles " << angles << " #cross.cor " << cc << "== " << get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)<< "\n"; + cout << "Allocation #Maps " << maps << " #Angles " << angles << " #cross.cor " << cc << "\n"; +//<< " == " << get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)<< "\n"; set_pointers(); } @@ -483,7 +512,8 @@ int bioem_RefMap::read_MRC(const char* filename,bioem_param& param) st2 = sqrt(st2 / float(nr * nc) - st * st); for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){ maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2; - //if(nmap+ntotRefMap==300) cout << i << " " << j << " " << nmap+ntotRefMap << " " << Ref[nmap+ntotRefMap].points[i][j] << "\n"; + //(nmap+ntotRefMap==300) + //cout <<"MAP:: " << i << " " << j << " " << maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] << "\n"; } } ntotRefMap += ns ; diff --git a/model.cpp b/model.cpp index de922c7af152dc538de95cea4acfbde4e45a25cd..3cb7733839f41bef2dd9de9d55d7241e5fc278dc 100644 --- a/model.cpp +++ b/model.cpp @@ -35,6 +35,8 @@ int bioem_model::readModel(const char* filemodel) ofstream exampleReadCoor; exampleReadCoor.open ("COORDREAD"); + + exampleReadCoor << "Text --- Number ---- x ---- y ---- z ---- radius ---- density\n"; int allocsize = 0; @@ -53,6 +55,17 @@ int bioem_model::readModel(const char* filemodel) int numres = 0; NormDen = 0.0; + string strfilename(filemodel); + + size_t foundpos= strfilename.find(".pdb"); + size_t endpos = strfilename.find_last_not_of(" \t"); + + if(foundpos > endpos){ + cout << "Warining:::: .pdb extension NOT dectected in file name \n"; + cout << "Warining:::: Are you sure you want to read a PDB? \n"; + } + + // cout << " HERE " << filemodel ; // for eachline in the file while (!input.eof()) @@ -127,7 +140,7 @@ int bioem_model::readModel(const char* filemodel) points[numres].point.pos[0] = (myfloat_t) x; points[numres].point.pos[1] = (myfloat_t) y; points[numres].point.pos[2] = (myfloat_t) z; - exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].density << "\n"; + exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].radius << " " << points[numres].density << "\n"; numres++; } } @@ -182,7 +195,7 @@ int bioem_model::readModel(const char* filemodel) points[numres].radius = (myfloat_t) tmpval[3]; points[numres].density = (myfloat_t) tmpval[4]; - exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].density << "\n"; + exampleReadCoor << "RESIDUE " << numres << " " << points[numres].point.pos[0] << " " << points[numres].point.pos[1] << " " << points[numres].point.pos[2] << " " << points[numres].radius << " " << points[numres].density << "\n"; NormDen += points[numres].density; numres++; } diff --git a/param.cpp b/param.cpp index 9f89b5977127851160d413a9c91f4f8f8680a442..16bcf317a7209137e8932537319b4fac661c7ada 100644 --- a/param.cpp +++ b/param.cpp @@ -77,7 +77,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) param_device.flipped=false; param_device.debugterm=false; param_device.writeCC=false; - param_device.CCwithBayes=false; + param_device.CCwithBayes=true; writeCTF=false; elecwavel=220; ignoreCCoff=false; @@ -236,16 +236,16 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) param_device.writeCC = true; cout << "Writing CrossCorrelations \n"; } - else if (strcmp(token, "CROSSCOR_DISPLACE") == 0) + else if (strcmp(token, "CROSSCOR_GRID_SPACE") == 0) { token = strtok(NULL, " "); param_device.CCdisplace=int(atoi(token)); if (param_device.CCdisplace < 0 ) { cout << "*** Error: Negative CROSSCOR_DISPLACE "; exit(1);} cout << "Writing Cross Correlation Displacement " << param_device.CCdisplace << " \n"; } - else if (strcmp(token, "CROSSCOR_WITHBAYESIAN") == 0) + else if (strcmp(token, "CROSSCOR_NOTBAYESIAN") == 0) { - param_device.CCwithBayes=true; + param_device.CCwithBayes=false; cout << "Using Bayesian Analysis to write Cross Correlation \n"; } else if (strcmp(token, "FLIPPED") == 0) //Key word if images are flipped for cross-correlation @@ -258,7 +258,7 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) ignoreCCoff = true; cout << "Ignoring Cross-Correlation offset \n"; } - else if (strcmp(token, "DO_AA_RADIUS") == 0) //If projecting CA with amino-acid radius + else if (strcmp(token, "PROJECT_RADIUS") == 0) //If projecting CA with amino-acid radius { doaaradius = true; cout << "Projecting corresponding radius \n"; @@ -300,8 +300,10 @@ int bioem_param::readParameters(const char* fileinput,const char* fileangles) // 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(!notuniformangles){ + 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);};