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
	      if (nummap % 128 == 0)
		{
		  cout << "..." << nummap << "\n";
		}
178
	      if(lasti+1 != param.param_device.NumberPixels && lastj+1 != 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
}