map.cpp 17.7 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

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

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

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
	numPixels = param.param_device.NumberPixels;
	refMapSize = param.param_device.NumberPixels * param.param_device.NumberPixels;
David Rohr's avatar
David Rohr committed
32
33
34
	// **************************************************************************************
	// ***********************Reading reference Particle Maps************************
	// **************************************************************************************
35
36
	int allocsize = 0;
	if (param.loadMap)
37
38
39
40
41
42
43
44
	{
		FILE* fp = fopen("maps.dump", "rb");
		if (fp == NULL)
		{
			cout << "Error opening dump file\n";
			exit(1);
		}
		fread(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
45
46
		maps = (myfloat_t*) mallocchk(ntotRefMap * refMapSize * sizeof(myfloat_t));
		fread(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
47
		fclose(fp);
48

David Rohr's avatar
David Rohr committed
49
		cout << "Particle Maps read from Map Dump\n";
50
	}
David Rohr's avatar
David Rohr committed
51
52
53
54
55
56
	else if(readMRC)
	{
		ntotRefMap=0;

		if(readMultMRC)
		{
57
			ifstream input(filemap);
David Rohr's avatar
David Rohr committed
58
59
60

			if (!input.good())
			{
61
				cout << "Failed to open file contaning MRC names: " << filemap << "\n";
David Rohr's avatar
David Rohr committed
62
63
64
65
66
67
68
69
70
71
				exit(0);
			}

			char line[512] = {0};
			char mapname[100];
			char tmpm[10] = {0};
			const char* indifile;

			while (!input.eof())
			{
David Rohr's avatar
David Rohr committed
72
73
				input.getline(line,511);
				char tmpVals[100]  = {0};
David Rohr's avatar
David Rohr committed
74

David Rohr's avatar
David Rohr committed
75
				strncpy (tmpVals,line,99);
David Rohr's avatar
David Rohr committed
76
77
78
				sscanf (tmpVals,"%99c",mapname);

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

			}
			cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
91
			cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
92
93
94
		}
		else
		{
95
			read_MRC(filemap,param);
David Rohr's avatar
David Rohr committed
96
			cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
97
			cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
David Rohr's avatar
David Rohr committed
98
99
		}
	}
100
101
	else
	{
102
103
104
		int nummap = -1;
		int lasti = 0;
		int lastj = 0;
105
		ifstream input(filemap);
106
107
108
109
110
111
		if (!input.good())
		{
			cout << "Particle Maps Failed to open file" << endl ;
			exit(1);
		}

David Rohr's avatar
David Rohr committed
112
113
		char line[512] = {0};
		char tmpLine[512] = {0};
114
                bool first=true;
115
116
117

		while (!input.eof())
		{
David Rohr's avatar
David Rohr committed
118
			input.getline(line, 511);
119

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

123
124
125
126
127
128
129
130
                        if(first){
			 	if (strcmp(token, "PARTICLE") != 0) {
                                 cout << "Missing correct Standard Map Format: PARTICLE HEADER\n"<< endl ;
	                        exit(1);
				}
                         first=false;
			}

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

David Rohr's avatar
David Rohr committed
159
				char tmpVals[36]  = {0};
160

161
162
				strncpy (tmpVals, line, 8);
				sscanf (tmpVals, "%d", &i);
163

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

167
168
				strncpy (tmpVals, line + 16, 12);
				sscanf (tmpVals, "%f", &z);
169
				//checking for Map limits
170
				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
177
178
179
180
181
182
183
				}
				else
				{
					cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
					exit(1);
				}
			}
		}
		cout << ".";
184
		ntotRefMap = nummap + 1;
185
		maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
David Rohr's avatar
David Rohr committed
186
		cout << "Particle Maps read from Standard File \n";
187

188
		if (param.dumpMap)
189
190
191
192
193
194
195
196
		{
			FILE* fp = fopen("maps.dump", "w+b");
			if (fp == NULL)
			{
				cout << "Error opening dump file\n";
				exit(1);
			}
			fwrite(&ntotRefMap, sizeof(ntotRefMap), 1, fp);
197
			fwrite(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
198
199
200
			fclose(fp);
		}
	}
David Rohr's avatar
David Rohr committed
201

202
203
204
205
	if (getenv("BIOEM_DEBUG_NMAPS"))
	{
	    ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS"));
	}
206

David Rohr's avatar
David Rohr committed
207
	cout << "Total Number of particles: " << ntotRefMap;
David Rohr's avatar
David Rohr committed
208
	cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
David Rohr's avatar
David Rohr committed
209

210
211
212
213
214
	return(0);
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
David Rohr's avatar
David Rohr committed
215
216
217
	// **************************************************************************************
	// ********** Routine that pre-calculates Kernels for Convolution **********************
	// ************************************************************************************
218

219
	RefMapsFFT = new mycomplex_t[ntotRefMap * param.FFTMapSize];
220

221
#pragma omp parallel for
222
223
	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
224
225
226
227
		const int num = omp_get_thread_num();
		myfloat_t* localMap = param.fft_scratch_real[num];
		mycomplex_t* localout = param.fft_scratch_complex[num];

228
		//Assigning localMap values to padded Map
229
		for(int i = 0; i < param.param_device.NumberPixels; i++)
230
		{
231
			for(int j = 0; j < param.param_device.NumberPixels; j++)
232
			{
233
				localMap[i * param.param_device.NumberPixels + j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
234
235
236
237
			}
		}

		//Calling FFT_Forward
238
		myfftw_execute_dft_r2c(param.fft_plan_r2c_forward, localMap, localout);
239

240
		//Saving the Reference CTFs (RefMap array possibly has incorrect alignment, so we copy here. Stupid but in fact does not matter.)
241
		mycomplex_t* RefMap = &RefMapsFFT[iRefMap * param.FFTMapSize];
242
		for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D ; i++ )
243
		{
244
245
			RefMap[i][0] = localout[i][0];
			RefMap[i][1] = localout[i][1];
246
247
248
249
		}
	}

	return(0);
250
}
251
252
253
254
255
256
257

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

258
259
260
	sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
	sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);

261
	//Precalculating cross-correlations of maps
262
	#pragma omp parallel for
263
264
	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
265
		myfloat_t sum, sumsquare;
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
		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);
}
282
283
284
285
286

void bioem_Probability::init(size_t maps, size_t angles, bioem& bio)
{
	nMaps = maps;
	nAngles = angles;
287
	ptr = bio.malloc_device_host(get_size(maps, angles, bio.param.param_device.writeAngles));
288
289
	set_pointers();
}
290

291
292
293
294
295
296
297
298
299
void bioem_Probability::initCC(size_t maps, size_t cc, bioem& bio)
{
        nMaps = maps;
        nCC = cc;
        ptr = bio.malloc_device_host(get_sizeCC(maps, cc, bio.param.param_device.writeCC));
        set_pointers();
}


David Rohr's avatar
David Rohr committed
300
301
302
303
304
305
void bioem_Probability::copyFrom(bioem_Probability* from, bioem& bio)
{
	bioem_Probability_map& pProbMap = getProbMap(0);
	bioem_Probability_map& pProbMapFrom = from->getProbMap(0);
	memcpy(&pProbMap, &pProbMapFrom, from->nMaps * sizeof(bioem_Probability_map));

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

David Rohr's avatar
David Rohr committed
317
318
int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
{
David Rohr's avatar
David Rohr committed
319
	myfloat_t st,st2;
David Rohr's avatar
David Rohr committed
320
321
322
323
324
325
326
	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;
327
	int  mx, my, mz,mapc, mapr, maps_local;
David Rohr's avatar
David Rohr committed
328
329
330
331
332
333
334
335
336
337
338
	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);

David Rohr's avatar
David Rohr committed
339
	if (n_range_viol0 < n_range_viol1) { //* guess endianism
David Rohr's avatar
David Rohr committed
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
		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);
370
	header_ok *= read_int(&maps_local,fin,swap);
David Rohr's avatar
David Rohr committed
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
	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);
398
	printf("     MAPS = %8d  (sections axis: 1=X,2=Y,3=Z)\n",maps_local);
David Rohr's avatar
David Rohr committed
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
	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);
	}

417
418
419
420
421
422
423
424
425
	if (ntotRefMap == 0)
	{
		maps = (myfloat_t*) mallocchk(refMapSize * sizeof(myfloat_t) * ns);
	}
	else
	{
		maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * (ntotRefMap + ns));
	}

David Rohr's avatar
David Rohr committed
426
427
428
429
430
	if(mode!=2)
	{
		cout << "ERROR with MRC mode " << mode << "\n";
		cout << "Currently mode 2 is the only one allowed" << "\n";
		exit(1);
431
432
	}
	else
David Rohr's avatar
David Rohr committed
433
	{
David Rohr's avatar
David Rohr committed
434
435
436
437
438
439
440
441
		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;
David Rohr's avatar
David Rohr committed
442
			exit(1);
David Rohr's avatar
David Rohr committed
443
444
445
446
447
448
449
450
451
		}

		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 ++ )
				{
452
453
					if (read_float(&currfloat,fin,swap)==0)
					{
David Rohr's avatar
David Rohr committed
454
455
						cout << "ERROR Converting Data: " <<  filename;
						exit(1);
456
457
458
459
					}
					else
					{
						maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = (myfloat_t) currfloat;
David Rohr's avatar
David Rohr committed
460
461
462
						st += currfloat;
						st2 += currfloat*currfloat;
					}
David Rohr's avatar
David Rohr committed
463
				}
David Rohr's avatar
David Rohr committed
464
465
466

				//Normaling maps to zero mean and unit standard deviation
				st /= float(nr*nc);
David Rohr's avatar
David Rohr committed
467
				st2 = sqrt(st2 / float(nr * nc) - st * st);
David Rohr's avatar
David Rohr committed
468
				for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
469
					maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2;
David Rohr's avatar
David Rohr committed
470
					//if(nmap+ntotRefMap==300) cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";
David Rohr's avatar
David Rohr committed
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
				}
		}
		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;
533
	int mapc, mapr, maps_local;
David Rohr's avatar
David Rohr committed
534
535
536
537
538
539
540
541
542
	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);
	}

David Rohr's avatar
David Rohr committed
543
	//* read header info
David Rohr's avatar
David Rohr committed
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
	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);
562
	header_ok *= read_int(&maps_local,fin,swap);
David Rohr's avatar
David Rohr committed
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
	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;
}