diff --git a/CMakeLists.txt b/CMakeLists.txt
index a97eed3837d6d84850f1c9b2697f4287e433d1c0..508c6fab630652cd4a6d4abcc2d67e64fbdd9899 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -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)
diff --git a/bioem.cpp b/bioem.cpp
index c5eeca600f63beb20ae1aa618e5e7a7637183347..70124a39d57552618b5d7c2432922d5f4c1bf768 100644
--- a/bioem.cpp
+++ b/bioem.cpp
@@ -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"))
{
diff --git a/include/map.h b/include/map.h
index f801d166810c446ab17e83ffb543c58be699abc1..5af0d9883aad073188362b172b7cc632e8881b4e 100644
--- a/include/map.h
+++ b/include/map.h
@@ -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]);}
diff --git a/main.cpp b/main.cpp
index cf8b518803b09de0c2ce9e5d8c0cc058cd34e471..587dac6ee7fcf2f4a64b6f01519c19c7ca54c5d7 100644
--- a/main.cpp
+++ b/main.cpp
@@ -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;
diff --git a/map.cpp b/map.cpp
index 7e1039751327323b1441c4c9aa531dd896f9f286..4ede0bd687392ab57d811deae913d4fe22ec4bf4 100644
--- a/map.cpp
+++ b/map.cpp
@@ -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];
+ 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 ae189ab9f12efb560e185feb97cb74f0192b469d..34e0ebf124c645c2ce0cdaa9b1b35db571c75e2d 100644
--- a/param.cpp
+++ b/param.cpp
@@ -37,9 +37,9 @@ bioem_param::bioem_param()
int bioem_param::readParameters()
{
- /**************************************************************************************/
- /***************************** Reading Input Parameters ******************************/
- /**************************************************************************************/
+ //**************************************************************************************
+ //***************************** Reading Input Parameters ******************************
+ //**************************************************************************************
ifstream input(fileinput);
if (!input.good())