map.cpp 18.2 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
				}
				else
				{
					cout << "PROBLEM READING MAP (Map number " << nummap << ", i " << i << ", j " << j << ")" << "\n";
					exit(1);
				}
			}
		}
Pilar Cossio's avatar
Pilar Cossio committed
183
184
185
186
187
                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);
                }
188
		cout << ".";
189
		ntotRefMap = nummap + 1;
190
		maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
David Rohr's avatar
David Rohr committed
191
		cout << "Particle Maps read from Standard File \n";
192

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

207
208
209
210
	if (getenv("BIOEM_DEBUG_NMAPS"))
	{
	    ntotRefMap = atoi(getenv("BIOEM_DEBUG_NMAPS"));
	}
211

David Rohr's avatar
David Rohr committed
212
	cout << "Total Number of particles: " << ntotRefMap;
David Rohr's avatar
David Rohr committed
213
	cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
David Rohr's avatar
David Rohr committed
214

215
216
217
218
219
	return(0);
}

int bioem_RefMap::PreCalculateMapsFFT(bioem_param& param)
{
David Rohr's avatar
David Rohr committed
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
	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
229
230
231
232
		const int num = omp_get_thread_num();
		myfloat_t* localMap = param.fft_scratch_real[num];
		mycomplex_t* localout = param.fft_scratch_complex[num];

233
		//Assigning localMap values to padded Map
234
		for(int i = 0; i < param.param_device.NumberPixels; i++)
235
		{
236
			for(int j = 0; j < param.param_device.NumberPixels; j++)
237
			{
238
				localMap[i * param.param_device.NumberPixels + j] = maps[iRefMap * refMapSize + i * param.param_device.NumberPixels + j];
239
240
241
242
			}
		}

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

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

	return(0);
255
}
256
257
258
259
260
261
262

int bioem_RefMap::precalculate(bioem_param& param, bioem& bio)
{
	// **************************************************************************************
	// *******************************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
	//Precalculating cross-correlations of maps
267
	#pragma omp parallel for
268
269
	for (int iRefMap = 0; iRefMap < ntotRefMap ; iRefMap++)
	{
270
		myfloat_t sum, sumsquare;
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
		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
void bioem_Probability::init(size_t maps, size_t angles, size_t cc, bioem& bio)
289
290
291
{
	nMaps = maps;
	nAngles = angles;
292
293
	nCC = cc;
	ptr = bio.malloc_device_host(get_size(maps, angles, cc, bio.param.param_device.writeAngles, bio.param.param_device.writeCC));
294
295
	set_pointers();
}
296

David Rohr's avatar
David Rohr committed
297
298
299
300
301
302
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));

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

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

David Rohr's avatar
David Rohr committed
324
325
int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
{
David Rohr's avatar
David Rohr committed
326
	myfloat_t st,st2;
David Rohr's avatar
David Rohr committed
327
328
329
330
331
332
333
	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;
334
	int  mx, my, mz,mapc, mapr, maps_local;
David Rohr's avatar
David Rohr committed
335
336
337
338
339
340
341
342
343
344
345
	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
346
	if (n_range_viol0 < n_range_viol1) { //* guess endianism
David Rohr's avatar
David Rohr committed
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
		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);
377
	header_ok *= read_int(&maps_local,fin,swap);
David Rohr's avatar
David Rohr committed
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
	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);
405
	printf("     MAPS = %8d  (sections axis: 1=X,2=Y,3=Z)\n",maps_local);
David Rohr's avatar
David Rohr committed
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
	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);
	}

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

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

				//Normaling maps to zero mean and unit standard deviation
				st /= float(nr*nc);
David Rohr's avatar
David Rohr committed
474
				st2 = sqrt(st2 / float(nr * nc) - st * st);
David Rohr's avatar
David Rohr committed
475
				for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
476
					maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] = maps[(nmap + ntotRefMap) * refMapSize + i * numPixels + j] / st2 - st/st2;
David Rohr's avatar
David Rohr committed
477
					//if(nmap+ntotRefMap==300) cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";
David Rohr's avatar
David Rohr committed
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
533
534
535
536
537
538
539
				}
		}
		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;
540
	int mapc, mapr, maps_local;
David Rohr's avatar
David Rohr committed
541
542
543
544
545
546
547
548
549
	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
550
	//* read header info
David Rohr's avatar
David Rohr committed
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
	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);
569
	header_ok *= read_int(&maps_local,fin,swap);
David Rohr's avatar
David Rohr committed
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
	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;
}