map.cpp 18.8 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
75
76
77
78
79
80
	  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};

	      strncpy (tmpVals,line,99);
	      sscanf (tmpVals,"%99c",mapname);

	      // Check for last line
	      strncpy (tmpm,mapname,3);
	      if(strcmp(tmpm,"XXX")!=0)
81
		{
82
83
84
		  indifile = mapname;
		  //Reading Multiple MRC
		  read_MRC(indifile,param);
85
		}
86
87
	      for(int i=0;i<3;i++)mapname[i] = 'X';
	      for(int i=3;i<100;i++)mapname[i] = 0;
88

89
90
91
92
93
94
95
96
97
	    }
	  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" ;
98
	}
99
100
101
102
103
104
105
106
    }
  else
    {
      int nummap = -1;
      int lasti = 0;
      int lastj = 0;
      ifstream input(filemap);
      if (!input.good())
David Rohr's avatar
David Rohr committed
107
	{
108
109
110
111
112
113
114
	  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
115

116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
      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;
	  }

	  if (strcmp(token, "PARTICLE") == 0) // to count the number of maps
	    {
	      nummap++;
	      if (allocsize == 0)
David Rohr's avatar
David Rohr committed
135
		{
136
137
		  allocsize = 64;
		  maps = (myfloat_t*) mallocchk(refMapSize * sizeof(myfloat_t) * allocsize);
138
		}
139
	      else if (nummap + 1 >= allocsize)
140
		{
141
142
		  allocsize *= 2;
		  maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * allocsize);
David Rohr's avatar
David Rohr committed
143
		}
144
145
146
147
148
	      if (nummap % 128 == 0)
		{
		  cout << "..." << nummap << "\n";
		}
	      if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels && nummap > 0)
149
		{
150
151
		  cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
		  exit(1);
152
		}
153
154
155
156
157
	    }
	  else
	    {
	      int i, j;
	      float z;
158

159
	      char tmpVals[36]  = {0};
160

161
162
163
164
165
166
167
168
169
170
	      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 < param.param_device.NumberPixels && j > 0 && j - 1 < param.param_device.NumberPixels)
171
		{
172
173
174
		  maps[nummap * refMapSize + (i - 1) * numPixels + (j - 1)] = (myfloat_t) z;
		  lasti = i;
		  lastj = j;
175
		}
176
	      else
177
		{
178
179
		  cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
		  exit(1);
180
		}
181
182
183
184
185
186
	    }
	}
      if(lasti != param.param_device.NumberPixels && lastj != param.param_device.NumberPixels )
	{
	  cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << lasti << ", j " << lastj << ")" << "\n";
	  exit(1);
187
	}
188
189
190
191
      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
192

193
      if (param.dumpMap)
194
	{
195
196
197
198
199
200
201
202
203
	  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);
204
	}
205
    }
206

207
208
209
210
  if (getenv("BIOEM_DEBUG_NMAPS"))
    {
      ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS"));
    }
David Rohr's avatar
David Rohr committed
211

212
213
214
215
  cout << "Total Number of particles: " << ntotRefMap;
  cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";

  return(0);
216
217
218
219
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
220
221
222
  // **************************************************************************************
  // ********** Routine that pre-calculates Kernels for Convolution **********************
  // ************************************************************************************
223

224
  RefMapsFFT = new mycomplex_t[ntotRefMap * param.FFTMapSize];
225

226
#pragma omp parallel for
227
228
229
230
231
232
233
234
  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++)
235
	{
236
237
238
239
240
	  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];
	    }
	}
241

242
243
      //Calling FFT_Forward
      myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
244

245
246
247
248
249
250
251
      //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]; 
252
	}
253
    }
254

255
  return(0);
256
}
257
258
259

int bioem_RefMap::precalculate(bioem_param& param, bioem& bio)
{
260
261
262
  // **************************************************************************************
  // *******************************Precalculating Routine for Maps************************
  // **************************************************************************************
263

264
265
  sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
  sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
266

267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
  //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);
287
}
288

289
void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio)
290
{
291
292
293
294
295
296
  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));
  cout << "INI " << maps << " " << angles << " " << cc << " " << get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC)<< "\n";
  set_pointers();
297
}
298

David Rohr's avatar
David Rohr committed
299
300
void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio)
{
301
302
303
  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
304

305
306
307
  if (bio.param.param_device.writeAngles)
    {
      for (int iOrient = 0; iOrient < nAngles; iOrient ++)
David Rohr's avatar
David Rohr committed
308
	{
309
310
311
	  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
312
	}
313
    }
314

315
316
317
  if (bio.param.param_device.writeCC)
    {
      for (int iCC = 0; iCC < nCC; iCC ++)
318
	{
319
320
321
	  bioem_Probability_cc& pProbCC = getProbCC(0, iCC);
	  bioem_Probability_cc& pProbCCFrom = from->getProbCC(0, iCC);
	  memcpy(&pProbCC, &pProbCCFrom, from->nMaps * sizeof(bioem_Probability_cc));
322
	}
323
    }
David Rohr's avatar
David Rohr committed
324
325
}

David Rohr's avatar
David Rohr committed
326
327
int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
{
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
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
  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
447
448
	}

449
450
451
      for (count=0; count<(unsigned long)nsymbt; ++count) if (read_char_float(&currfloat,fin)==0) {
	  cout << "ERROR Converting Data: " <<  filename;
	  exit(1);
452
453
	}

454
      for ( int nmap = 0 ; nmap < ns ; nmap ++ )
David Rohr's avatar
David Rohr committed
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
	  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;
	      //if(nmap+ntotRefMap==300) cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";
	    }
David Rohr's avatar
David Rohr committed
481
	}
482
483
484
485
      ntotRefMap += ns ;
      //  cout << ntotRefMap << "\n";
    }
  fclose (fin);
David Rohr's avatar
David Rohr committed
486

487
  return(0);
David Rohr's avatar
David Rohr committed
488
489
490
}

int bioem_RefMap::read_float(float *currfloat, FILE *fin, int swap) {
491
492
493
494
495
496
497
498
499
500
501
502
503
  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
504
505
506
}

int  bioem_RefMap::read_int(int *currlong, FILE *fin, int swap) {
507
508
509
510
511
512
513
514
515
516
517
518
519
  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
520
521
}
int  bioem_RefMap::read_float_empty (FILE *fin) {
522
  float currfloat;
David Rohr's avatar
David Rohr committed
523

524
525
  if (fread(&currfloat,4,1,fin)!=1) return 0;
  return 1;
David Rohr's avatar
David Rohr committed
526
527
528
}

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

531
532
533
  if (fread(&currchar,1,1,fin)!=1) return 0;
  *currfloat=(float)currchar;
  return 1;
David Rohr's avatar
David Rohr committed
534
535
536
}

int  bioem_RefMap::test_mrc (const char *vol_file, int swap) {
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
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
  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
599
}