map.cpp 20 KB
Newer Older
Pilar Cossio's avatar
License  
Pilar Cossio committed
1
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 3 4
   < BioEM software for Bayesian inference of Electron Microscopy images>
   Copyright (C) 2014 Pilar Cossio, David Rohr and Gerhard Hummer.
   Max Planck Institute of Biophysics, Frankfurt, Germany.
David Rohr's avatar
David Rohr committed
5

6
   See license statement for terms of distribution.
Pilar Cossio's avatar
License  
Pilar Cossio committed
7 8 9

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

10 11 12 13 14
#include <fstream>
#include <iostream>
#include <stdio.h>
#include <stdlib.h>
#include <cstring>
15 16
#include <math.h>
#include <fftw3.h>
17

18 19 20 21
#ifdef WITH_OPENMP
#include <omp.h>
#endif

22
#include "map.h"
23
#include "param.h"
24
#include "bioem.h"
25

26 27
using namespace std;

28
int bioem_RefMap::readRefMaps(bioem_param& param, const char* filemap)
29
{
30 31 32 33 34 35 36 37 38 39
  numPixels = param.param_device.NumberPixels;
  refMapSize = param.param_device.NumberPixels * param.param_device.NumberPixels;
  // **************************************************************************************
  // ***********************Reading reference Particle Maps************************
  // **************************************************************************************
  int allocsize = 0;
  if (param.loadMap)
    {
      FILE* fp = fopen("maps.dump", "rb");
      if (fp == NULL)
40
	{
41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
	  cout << "Error opening dump file\n";
	  exit(1);
	}
      fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
      maps = (myfloat_t*) mallocchk(ntotRefMap * refMapSize * sizeof(myfloat_t));
      fread(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
      fclose(fp);

      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,511);
	      char tmpVals[100]  = {0};

75 76 77 78 79
		string strline(line);

 	//	 cout << "MRC File name:" << strline << "\n";


80 81 82 83 84
	      strncpy (tmpVals,line,99);
	      sscanf (tmpVals,"%99c",mapname);

	      // Check for last line
	      strncpy (tmpm,mapname,3);
85

86
	      if(strcmp(tmpm,"XXX")!=0)
87
		{
88
		indifile=strline.c_str();
Pilar Cossio's avatar
Pilar Cossio committed
89 90 91 92 93 94 95 96 97
          
                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";
                }

98 99
		  //Reading Multiple MRC
		  read_MRC(indifile,param);
100
		}
101 102
	      for(int i=0;i<3;i++)mapname[i] = 'X';
	      for(int i=3;i<100;i++)mapname[i] = 0;
103

104 105 106 107 108 109
	    }
	  cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
	  cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
	}
      else
	{
Pilar Cossio's avatar
Pilar Cossio committed
110 111 112 113 114 115 116 117 118 119 120

		 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";
                }

121 122 123
	  read_MRC(filemap,param);
	  cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
	  cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
124
	}
125 126 127 128 129 130 131 132
    }
  else
    {
      int nummap = -1;
      int lasti = 0;
      int lastj = 0;
      ifstream input(filemap);
      if (!input.good())
David Rohr's avatar
David Rohr committed
133
	{
134 135 136 137 138 139 140
	  cout << "Particle Maps Failed to open file" << endl ;
	  exit(1);
	}

      char line[512] = {0};
      char tmpLine[512] = {0};
      bool first=true;
David Rohr's avatar
David Rohr committed
141

Pilar Cossio's avatar
Pilar Cossio committed
142 143
      int countpix;

144 145 146 147 148 149 150 151 152 153 154 155 156 157 158
      while (!input.eof())
	{
	  input.getline(line, 511);

	  strncpy(tmpLine, line, strlen(line));
	  char *token = strtok(tmpLine, " ");

	  if(first){
	    if (strcmp(token, "PARTICLE") != 0) {
	      cout << "Missing correct Standard Map Format: PARTICLE HEADER\n"<< endl ;
	      exit(1);
	    }
	    first=false;
	  }

Pilar Cossio's avatar
Pilar Cossio committed
159

160 161 162
	  if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
	    {
	      nummap++;
Pilar Cossio's avatar
Pilar Cossio committed
163
	      countpix=0;
164
	      if (allocsize == 0)
David Rohr's avatar
David Rohr committed
165
		{
166 167
		  allocsize = 64;
		  maps = (myfloat_t*) mallocchk(refMapSize * sizeof(myfloat_t) * allocsize);
168
		}
169
	      else if (nummap + 1 >= allocsize)
170
		{
171 172
		  allocsize *= 2;
		  maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * allocsize);
David Rohr's avatar
David Rohr committed
173
		}
174 175 176 177 178
	      if (nummap % 128 == 0)
		{
		  cout << "..." << nummap << "\n";
		}
	      if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
179
		{
180 181
		  cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
		  exit(1);
182
		}
183 184 185 186 187
	    }
	  else
	    {
	      int i, j;
	      float z;
Pilar Cossio's avatar
Pilar Cossio committed
188
	      
189
	      char tmpVals[36]  = {0};
190

191 192 193 194 195 196
	      strncpy (tmpVals, line, 8);
	      sscanf (tmpVals, "%d", &i);

	      strncpy (tmpVals, line + 8, 8);
	      sscanf (tmpVals, "%d", &j);

Pilar Cossio's avatar
Pilar Cossio committed
197
	      strncpy (tmpVals, line + 16, 16);
198 199
	      sscanf (tmpVals, "%f", &z);
	      //checking for Map limits
Pilar Cossio's avatar
Pilar Cossio committed
200
	      if(i > -1 && i  < param.param_device.NumberPixels && j > -1 && j  < param.param_device.NumberPixels)
201
		{
Pilar Cossio's avatar
Pilar Cossio committed
202 203
		  countpix++;
		  maps[nummap * refMapSize + i * numPixels + j] = (myfloat_t) z;
204 205
		  lasti = i;
		  lastj = j;
Pilar Cossio's avatar
Pilar Cossio committed
206
	 //	 cout << countpix << " " << param.param_device.NumberPixels*param.param_device.NumberPixels << "\n";
207
		}
208
	      else
209
		{
210 211
		  cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
		  exit(1);
212
		}
213 214
	    }
	}
Pilar Cossio's avatar
Pilar Cossio committed
215
      if(lasti != param.param_device.NumberPixels-1 || lastj != param.param_device.NumberPixels-1 || countpix != param.param_device.NumberPixels*param.param_device.NumberPixels +1 )
216 217 218
	{
	  cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
	  exit(1);
219
	}
220 221 222 223
      cout << ".";
      ntotRefMap = nummap + 1;
      maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
      cout << "Particle Maps read from Standard File \n";
David Rohr's avatar
David Rohr committed
224

225
    }
226

Pilar Cossio's avatar
Pilar Cossio committed
227 228 229 230 231 232 233 234 235 236 237 238 239 240 241
//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);
        }


242 243 244 245
  if (getenv("BIOEM_DEBUG_NMAPS"))
    {
      ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS"));
    }
David Rohr's avatar
David Rohr committed
246

247 248 249 250
  cout << "Total Number of particles: " << ntotRefMap;
  cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";

  return(0);
251 252 253 254
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
255 256 257
  // **************************************************************************************
  // ********** Routine that pre-calculates Kernels for Convolution **********************
  // ************************************************************************************
258

259
  RefMapsFFT = new mycomplex_t[ntotRefMap * param.FFTMapSize];
260

261
#pragma omp parallel for
262 263 264 265 266 267 268 269
  for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
    {
      const int num = omp_get_thread_num();
      myfloat_t* localMap = param.fft_scratch_real[num];
      mycomplex_t* localout = param.fft_scratch_complex[num];

      //Assigning localMap values to padded Map
      for(int i = 0; i < param.param_device.NumberPixels; i++)
270
	{
271 272 273 274 275
	  for(int j = 0; j < param.param_device.NumberPixels; j++)
	    {
	      localMap[i * param.param_device.NumberPixels + j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
	    }
	}
276

277 278
      //Calling FFT_Forward
      myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
279

280 281 282 283 284 285 286
      //Saving the Reference CTFs (RefMap array possibly has incorrect alignment, so we copy here. Stupid but in fact does not matter.)
      mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.FFTMapSize];
               
      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]; 
287
	}
288
    }
289

290
  return(0);
291
}
292 293 294

int bioem_RefMap::precalculate(bioem_param& param, bioem& bio)
{
295 296 297
  // **************************************************************************************
  // *******************************Precalculating Routine for Maps************************
  // **************************************************************************************
298

299 300
  sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
  sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
301

302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321
  //Precalculating cross-correlations of maps
#pragma omp parallel for
  for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
    {
      myfloat_t sum, sumsquare;
      bio.calcross_cor(getmap(iRefMap), sum, sumsquare);
      //Storing Crosscorrelations in Map class
      sum_RefMap[iRefMap] = sum;
      sumsquare_RefMap[iRefMap] = sumsquare;
    }

  // Precalculating Maps in Fourier space
  if (bio.FFTAlgo)
    {
      PreCalculateMapsFFT(param);
      free(maps);
      maps = NULL;
    }

  return(0);
322
}
323

324
void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio)
325
{
326 327 328 329
  nMaps = maps;
  nAngles = angles;
  nCC = cc;
  ptr = bio.malloc_device_host(get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC));
Pilar Cossio's avatar
Pilar Cossio committed
330 331
  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";
332
  set_pointers();
333
}
334

David Rohr's avatar
David Rohr committed
335 336
void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio)
{
337 338 339
  bioem_Probability_map& pProbMap = getProbMap(0);
  bioem_Probability_map& pProbMapFrom = from->getProbMap(0);
  memcpy(&pProbMap, &pProbMapFrom, from->nMaps * sizeof(bioem_Probability_map));
David Rohr's avatar
David Rohr committed
340

341 342 343
  if (bio.param.param_device.writeAngles)
    {
      for (int iOrient = 0; iOrient < nAngles; iOrient ++)
David Rohr's avatar
David Rohr committed
344
	{
345 346 347
	  bioem_Probability_angle& pProbAngle = getProbAngle(0, iOrient);
	  bioem_Probability_angle& pProbAngleFrom = from->getProbAngle(0, iOrient);
	  memcpy(&pProbAngle, &pProbAngleFrom, from->nMaps * sizeof(bioem_Probability_angle));
David Rohr's avatar
David Rohr committed
348
	}
349
    }
350

351 352 353
  if (bio.param.param_device.writeCC)
    {
      for (int iCC = 0; iCC < nCC; iCC ++)
354
	{
355 356 357
	  bioem_Probability_cc& pProbCC = getProbCC(0, iCC);
	  bioem_Probability_cc& pProbCCFrom = from->getProbCC(0, iCC);
	  memcpy(&pProbCC, &pProbCCFrom, from->nMaps * sizeof(bioem_Probability_cc));
358
	}
359
    }
David Rohr's avatar
David Rohr committed
360 361
}

David Rohr's avatar
David Rohr committed
362 363
int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
{
364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482
  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_local;
  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_local,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_local);
     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 (ntotRefMap == 0)
    {
      maps = (myfloat_t*) mallocchk(refMapSize * sizeof(myfloat_t) * ns);
    }
  else
    {
      maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * (ntotRefMap + ns));
    }

  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);
David Rohr's avatar
David Rohr committed
483 484
	}

485 486 487
      for (count=0; count<(unsigned long)nsymbt; ++count) if (read_char_float(&currfloat,fin)==0) {
	  cout << "ERROR Converting Data: " <<  filename;
	  exit(1);
488 489
	}

490
      for ( int nmap = 0 ; nmap < ns ; nmap ++ )
David Rohr's avatar
David Rohr committed
491
	{
492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514
	  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
		  {
		    maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = (myfloat_t) 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 ++ ){
	      maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2;
Pilar Cossio's avatar
Pilar Cossio committed
515 516
	      //(nmap+ntotRefMap==300) 
	      //cout <<"MAP:: " << i << " " << j << " " <<  maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j]  << "\n";
517
	    }
David Rohr's avatar
David Rohr committed
518
	}
519 520 521 522
      ntotRefMap += ns ;
      //  cout << ntotRefMap << "\n";
    }
  fclose (fin);
David Rohr's avatar
David Rohr committed
523

524
  return(0);
David Rohr's avatar
David Rohr committed
525 526 527
}

int bioem_RefMap::read_float(float *currfloat, FILE *fin, int swap) {
528 529 530 531 532 533 534 535 536 537 538 539 540
  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;
David Rohr's avatar
David Rohr committed
541 542 543
}

int  bioem_RefMap::read_int(int *currlong, FILE *fin, int swap) {
544 545 546 547 548 549 550 551 552 553 554 555 556
  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;
David Rohr's avatar
David Rohr committed
557 558
}
int  bioem_RefMap::read_float_empty (FILE *fin) {
559
  float currfloat;
David Rohr's avatar
David Rohr committed
560

561 562
  if (fread(&currfloat,4,1,fin)!=1) return 0;
  return 1;
David Rohr's avatar
David Rohr committed
563 564 565
}

int  bioem_RefMap::read_char_float (float *currfloat, FILE *fin) {
566
  char currchar;
David Rohr's avatar
David Rohr committed
567

568 569 570
  if (fread(&currchar,1,1,fin)!=1) return 0;
  *currfloat=(float)currchar;
  return 1;
David Rohr's avatar
David Rohr committed
571 572 573
}

int  bioem_RefMap::test_mrc (const char *vol_file, int swap) {
574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635
  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_local;
  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_local,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;
David Rohr's avatar
David Rohr committed
636
}