EMAN2
mrcio.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#include <cstring>
33#include <climits>
34
35#include <sys/stat.h>
36
37#include "mrcio.h"
38#include "portable_fileio.h"
39#include "geometry.h"
40#include "util.h"
41#include "ctf.h"
42#include "transform.h"
43
44using namespace EMAN;
45
46const char *MrcIO::CTF_MAGIC = "!-";
47const char *MrcIO::SHORT_CTF_MAGIC = "!$";
48
49MrcIO::MrcIO(const string & fname, IOMode rw)
50: ImageIO(fname, rw), mode_size(0),
51 isFEI(false), is_ri(0), is_new_file(false),
52 is_transpose(false), is_stack(false), stack_size(1),
53 is_8_bit_packed(false), use_given_dimensions(true),
54 rendermin(0.0), rendermax(0.0), renderbits(16)
55{
56 memset(&mrch, 0, sizeof(MrcHeader));
57 is_big_endian = ByteOrder::is_host_big_endian();
58}
59
60MrcIO::~MrcIO()
61{
62 if (file) {
63 fclose(file);
64 file = NULL;
65 }
66}
67
68void MrcIO::init()
69{
70 ENTERFUNC;
71
72 if (initialized) {
73 return;
74 }
75
76 setbuf (stdout, NULL);
77
78 IOMode rwmode;
79
80 if (rw_mode == WRITE_ONLY) {
81 rwmode = READ_WRITE;
82 }
83 else {
84 rwmode = rw_mode;
85 }
86
87 int error_type;
88 struct stat status;
89
90 error_type = stat(filename.c_str(), & status);
91 is_new_file = (error_type != 0);
92
93 initialized = true;
94
95// mrcfile = sfopen(filename, rwmode, &is_new_file);
96 file = sfopen(filename, rwmode, NULL);
97
98 string ext = Util::get_filename_ext(filename);
99
100 if (ext != "") {
101 if (ext == "raw" || ext == "RAW") {
102 isFEI = true;
103 }
104
105 if (ext == "mrcs" || ext == "MRCS") {
106 is_stack = true;
107 }
108 }
109
110 if (! is_new_file) {
111 if (fread(&mrch, sizeof(MrcHeader), 1, file) != 1) {
112 throw ImageReadException(filename, "MRC header");
113 }
114
115 bool do_swap, have_err;
116
117 check_swap((const int *) (& mrch), filename.c_str(), true,
118 do_swap, have_err);
119
120 if (have_err || ! is_valid(&mrch)) {
121 throw ImageReadException(filename, "invalid MRC");
122 }
123
124 if (do_swap) {
125 swap_header(mrch);
126
127 is_big_endian = ! ByteOrder::is_host_big_endian();
128 }
129 else {
130 is_big_endian = ByteOrder::is_host_big_endian();
131 }
132
133 int max_labels = Util::get_min(mrch.nlabels, (int) MRC_NUM_LABELS);
134
135 for (int ilabel = 0; ilabel < max_labels; ilabel++) {
136 Util::replace_non_ascii(mrch.labels[ilabel], MRC_LABEL_SIZE);
137 }
138
139 // become_host_endian((int *) &mrch, NUM_4BYTES_PRE_MAP);
140 // become_host_endian((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
141 mode_size = get_mode_size(mrch.mode);
142
143 if (is_complex_mode()) {
144 is_ri = 1;
145 }
146
147 if (mrch.nxstart != 0 || mrch.nystart != 0 || mrch.nzstart != 0) {
148 LOGWARN("nx/ny/nz start not zero");
149 }
150
151 if (is_complex_mode()) {
152 mrch.nx *= 2;
153 }
154
155 if (mrch.xlen == 0) {
156 mrch.xlen = 1.0;
157 }
158
159 if (mrch.ylen == 0) {
160 mrch.ylen = 1.0;
161 }
162
163 if (mrch.zlen == 0) {
164 mrch.zlen = 1.0;
165 }
166
167 if (mrch.nlabels > 0) {
168 if (string(mrch.labels[0],3) == "Fei") {
169 isFEI = true;
170 }
171 }
172
173 if (mrch.mapc == 2 && mrch.mapr == 1) {
174 is_transpose = true;
175 }
176
177 if (is_stack) {
178 stack_size = mrch.nz;
179 mrch.nz = 1;
180 }
181
182 // Stuff added for 8 bit packed mode,
183 // with 2 4-bit values packed into each 8-bit byte:
184 float ny_to_nx_ratio;
185
186 if (mrch.nx > 0) {
187 ny_to_nx_ratio = (float) mrch.ny / (float) mrch.nx;
188 }
189 else {
190 ny_to_nx_ratio = 1.0;
191 }
192
193 bool have_packed_label = (mrch.nlabels > 0 &&
194 strstr(mrch.labels[0], "4 bits packed") != NULL);
195 bool have_packed_filename =
196 (strstr(filename.c_str(), "4_bits_packed") != NULL);
197
198 bool ny_twice_nx = (fabs(ny_to_nx_ratio - 2.0f) <= 0.1f);
199 bool double_nx = (ny_twice_nx && mrch.mode == MRC_UCHAR &&
200 (have_packed_label || have_packed_filename));
201
202 use_given_dimensions = (! double_nx);
203
204 is_8_bit_packed = (mrch.mode == MRC_UHEX || (mrch.imod_flags & 16) || double_nx);
205
206 if (getenv("DISALLOW_PACKED_FORMAT")) {
207 use_given_dimensions = true;
208 is_8_bit_packed = false;
209 }
210
211 if (is_8_bit_packed) {
212 mrch.mode = MRC_UHEX;
213
214 if (! use_given_dimensions) {
215 mrch.nx = mrch.nx * 2;
216 }
217 }
218 }
219
220 EXITFUNC;
221}
222
223bool MrcIO::is_image_big_endian()
224{
225 init();
226
227 return is_big_endian;
228}
229
230void MrcIO::check_swap(const int * data, const char * filnam, bool show_errors,
231 bool & do_swap, bool & have_err)
232{
233 int nx = data[0];
234 int ny = data[1];
235 int nz = data[2];
236 int mrcmode = data[3];
237 int mapc = data[16]; // Which axis corresponds to Columns (X=1, Y=2, Z=3)
238 int mapr = data[17]; // Which axis corresponds to Rows (X=1, Y=2, Z=3)
239 int maps = data[18]; // Which axis corresponds to Sections (X=1, Y=2, Z=3)
240 int nsymbt = data[23]; // the extra bytes for symmetry information
241 int mach = data[44]; // Machine endian stamp
242
243 int nxw, nyw, nzw, modew, mapcw, maprw, mapsw, nsymw, machw;
244
245 nxw = nx;
246 nyw = ny;
247 nzw = nz;
248 modew = mrcmode;
249 mapcw = mapc;
250 maprw = mapr;
251 mapsw = maps;
252 nsymw = nsymbt;
253 machw = mach;
254
255 ByteOrder::swap_bytes(&nxw);
256 ByteOrder::swap_bytes(&nyw);
257 ByteOrder::swap_bytes(&nzw);
258 ByteOrder::swap_bytes(&modew);
259 ByteOrder::swap_bytes(&mapcw);
260 ByteOrder::swap_bytes(&maprw);
261 ByteOrder::swap_bytes(&mapsw);
262 ByteOrder::swap_bytes(&nsymw);
263 ByteOrder::swap_bytes(&machw);
264
265 const int max_dim = 1 << 20;
266
267 int actual_stamp = MrcIO::generate_machine_stamp();
268
269 bool debug = (getenv("DEBUG_MRC_SANITY"));
270
271 string errmsg;
272
273 have_err = false;
274
275 if (mach == actual_stamp) {
276 do_swap = false;
277 }
278 else if (machw == actual_stamp) {
279 do_swap = true;
280 }
281 else {
282 if (mrcmode == 0) {
283 if (1 <= mapr && mapr <= 3 &&
284 1 <= mapc && mapc <= 3 &&
285 1 <= maps && maps <= 3) {
286 do_swap = false;
287 }
288 else if (1 <= maprw && maprw <= 3 &&
289 1 <= mapcw && mapcw <= 3 &&
290 1 <= mapsw && mapsw <= 3) {
291 do_swap = true;
292 }
293 else {
294 double ave_xyz = ((double) nx + (double) ny + (double) nz) / 3.0;
295 double ave_xyzw = ((double) nxw + (double) nyw + (double) nzw) / 3.0;
296
297 if (nx > 0 && ny > 0 && nz > 0 && ave_xyz <= max_dim) {
298 do_swap = false;
299 }
300 else if (nxw > 0 && nyw > 0 && nzw > 0 && ave_xyzw <= max_dim) {
301 do_swap = true;
302 }
303 else {
304 have_err = true;
305 do_swap = false;
306 errmsg = "MRC image dimensions nonpositive or too large.";
307 }
308 }
309 }
310 else {
311 if (mrcmode > 0 && mrcmode < 128) {
312 do_swap = false;
313 }
314 else if (modew > 0 && modew < 128) {
315 do_swap = true;
316 }
317 else {
318 have_err = true;
319 do_swap = false;
320 errmsg = "MRC mode is not from 0 to 127.";
321 }
322 }
323 }
324
325 if (debug || (have_err && show_errors)) {
326 if (filnam) {
327 printf ("file name = '%s'.\n", filnam);
328 }
329
330 printf ("stamp: mach, file, swapd = %8.0x %8.0x %8.0x\n",
331 actual_stamp, mach, machw);
332 printf ("stamp: mach, file, swapd = %d %d %d\n",
333 actual_stamp, mach, machw);
334 printf ("nx, ny, nz, mode = %d %d %d %d\n", nx, ny, nz, mrcmode);
335 printf ("nx, ny, nz, mode swapped = %d %d %d %d\n", nxw, nyw, nzw, modew);
336 printf ("mapc, mapr, maps = %d %d %d\n", mapc, mapr, maps);
337 printf ("mapc, mapr, maps swapped = %d %d %d\n", mapcw, maprw, mapsw);
338 printf ("%s\n", errmsg.c_str());
339 }
340}
341
342bool MrcIO::is_valid(const void * first_block, off_t file_size)
343{
344 ENTERFUNC;
345
346 if (! first_block) { // Null pointer - no data
347 return false;
348 }
349
350 const int * data = (const int *) first_block;
351
352 int nx = data[0];
353 int ny = data[1];
354 int nz = data[2];
355 int mrcmode = data[3];
356 int nsymbt = data[23]; // the extra bytes for symmetry information
357
358 const int max_dim = 1 << 20;
359
360 bool do_swap, have_err;
361
362 check_swap(data, NULL, false, do_swap, have_err);
363
364 if (have_err) {
365 return false;
366 }
367
368 if (do_swap) {
369 ByteOrder::swap_bytes(&nx);
370 ByteOrder::swap_bytes(&ny);
371 ByteOrder::swap_bytes(&nz);
372 ByteOrder::swap_bytes(&mrcmode);
373 ByteOrder::swap_bytes(&nsymbt);
374 }
375
376 if (mrcmode == MRC_SHORT_COMPLEX || mrcmode == MRC_FLOAT_COMPLEX) {
377 nx *= 2;
378 }
379
380 if ((mrcmode >= MRC_UCHAR &&
381 (mrcmode < MRC_UNKNOWN || mrcmode == MRC_UHEX)) &&
382 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim)
383 && (nz > 0 && nz < max_dim)) {
384
385//#ifndef SPIDERMRC // Spider MRC files don't satisfy the following test
386
387 if (file_size > 0) {
388 off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz *
389 (off_t)get_mode_size(mrcmode) +
390 (off_t)sizeof(MrcHeader) + nsymbt;
391
392 if (file_size == file_size1) {
393 return true;
394 }
395
396// return false;
397
398 // when size doesn't match, print error message instead of make it fail
399
400 LOGWARN("image size check fails, still try to read it...");
401 }
402 else {
403 return true;
404 }
405
406//#endif // SPIDERMRC
407
408 return true;
409 }
410
411 EXITFUNC;
412
413 return false;
414}
415
416int MrcIO::read_header(Dict & dict, int image_index, const Region * area, bool is_3d)
417{
418 init();
419
420 if (isFEI) {
421 return read_fei_header(dict, image_index, area, is_3d);
422 }
423 else {
424 return read_mrc_header(dict, image_index, area, is_3d);
425 }
426}
427
428int MrcIO::read_mrc_header(Dict & dict, int image_index, const Region * area, bool)
429{
430 ENTERFUNC;
431
432 if (image_index < 0) {
433 image_index = 0;
434 }
435
436 if (image_index != 0 && ! is_stack) {
437 throw ImageReadException(filename,
438 "no stack allowed for MRC image. For take 2D slice out of 3D image, "
439 "read the 3D image first, then use get_clip().");
440 }
441
442 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
443
444 int xlen = 0, ylen = 0, zlen = 0;
445
446 EMUtil::get_region_dims(area, mrch.nx, & xlen, mrch.ny, & ylen, mrch.nz, & zlen);
447
448 dict["nx"] = xlen;
449 dict["ny"] = ylen;
450 dict["nz"] = zlen;
451 dict["MRC.nx"] = mrch.nx;
452 dict["MRC.ny"] = mrch.ny;
453 dict["MRC.nz"] = mrch.nz;
454
455 dict["datatype"] = to_em_datatype(mrch.mode);
456
457 dict["MRC.nxstart"] = mrch.nxstart;
458 dict["MRC.nystart"] = mrch.nystart;
459 dict["MRC.nzstart"] = mrch.nzstart;
460
461 dict["MRC.mx"] = mrch.mx;
462 dict["MRC.my"] = mrch.my;
463 dict["MRC.mz"] = mrch.mz;
464
465 dict["MRC.xlen"] = mrch.xlen;
466 dict["MRC.ylen"] = mrch.ylen;
467 dict["MRC.zlen"] = mrch.zlen;
468
469 dict["MRC.alpha"] = mrch.alpha;
470 dict["MRC.beta"] = mrch.beta;
471 dict["MRC.gamma"] = mrch.gamma;
472
473 dict["MRC.mapc"] = mrch.mapc;
474 dict["MRC.mapr"] = mrch.mapr;
475 dict["MRC.maps"] = mrch.maps;
476
477 dict["MRC.minimum"] = mrch.amin;
478 dict["MRC.maximum"] = mrch.amax;
479 dict["MRC.mean"] = mrch.amean;
480// dict["minimum"] = mrch.amin;
481// dict["maximum"] = mrch.amax;
482// dict["mean"] = mrch.amean;
483
484 dict["MRC.ispg"] = mrch.ispg;
485 dict["MRC.nsymbt"] = mrch.nsymbt;
486
487 float apx = mrch.xlen / mrch.mx;
488 float apy = mrch.ylen / mrch.my;
489 float apz = mrch.zlen / mrch.mz;
490
491 if (apx > 1000 || apx < 0.01) {
492 dict["apix_x"] = 1.0f;
493 }
494 else {
495 dict["apix_x"] = apx;
496 }
497
498 if (apy > 1000 || apy < 0.01) {
499 dict["apix_y"] = 1.0f;
500 }
501 else {
502 dict["apix_y"] = apy;
503 }
504
505 if (apz > 1000 || apz < 0.01) {
506 dict["apix_z"] = 1.0f;
507 }
508 else {
509 dict["apix_z"] = apz;
510 }
511
512 if (area) {
513 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
514 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
515
516 if (area->get_ndim() == 3 && mrch.nz > 1) {
517 dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
518 }
519 else {
520 dict["origin_z"] = mrch.zorigin;
521 }
522 }
523 else {
524 dict["origin_x"] = mrch.xorigin;
525 dict["origin_y"] = mrch.yorigin;
526 dict["origin_z"] = mrch.zorigin;
527 }
528
529 if (is_complex_mode()) {
530 dict["is_complex"] = 1;
531 dict["is_complex_ri"] = 1;
532 }
533
534 dict["MRC.machinestamp"] = mrch.machinestamp;
535
536 dict["MRC.rms"] = mrch.rms;
537 dict["sigma"] = mrch.rms;
538 dict["MRC.nlabels"] = mrch.nlabels;
539
540 char label[32];
541
542 for (int i = 0; i < mrch.nlabels; i++) {
543 sprintf(label, "MRC.label%d", i);
544 dict[string(label)] = string(mrch.labels[i], MRC_LABEL_SIZE);
545 }
546
547 EMAN1Ctf ctf_;
548
549 if (read_ctf(ctf_) == 0) {
550 vector<float> vctf = ctf_.to_vector();
551 dict["ctf"] = vctf;
552 }
553
554 if (is_transpose) {
555 dict["nx"] = ylen;
556 dict["ny"] = xlen;
557 dict["MRC.nx"] = mrch.ny;
558 dict["MRC.ny"] = mrch.nx;
559 dict["MRC.mx"] = mrch.my;
560 dict["MRC.my"] = mrch.mx;
561 dict["apix_x"] = mrch.ylen / mrch.my;
562 dict["apix_y"] = mrch.xlen / mrch.mx;
563 dict["origin_x"] = mrch.yorigin;
564 dict["origin_y"] = mrch.xorigin;
565 dict["MRC.nxstart"] = mrch.nystart;
566 dict["MRC.nystart"] = mrch.nxstart;
567 }
568
569 EXITFUNC;
570
571 return 0;
572}
573
574int MrcIO::read_fei_header(Dict & dict, int image_index, const Region * area, bool)
575{
576 ENTERFUNC;
577
578 if(image_index < 0) {
579 image_index = 0;
580 }
581
582 init();
583
584 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
585
586 int xlen = 0, ylen = 0, zlen = 0;
587
588 EMUtil::get_region_dims(area, feimrch.nx, & xlen, feimrch.ny, & ylen, feimrch.nz, & zlen);
589
590 dict["nx"] = xlen;
591 dict["ny"] = ylen;
592 dict["nz"] = zlen;
593 dict["FEIMRC.nx"] = feimrch.nx;
594 dict["FEIMRC.ny"] = feimrch.ny;
595 dict["FEIMRC.nz"] = feimrch.nz;
596
597 dict["datatype"] = to_em_datatype(feimrch.mode); //=1, FEI-MRC file always use short for data type
598
599 dict["FEIMRC.nxstart"] = feimrch.nxstart;
600 dict["FEIMRC.nystart"] = feimrch.nystart;
601 dict["FEIMRC.nzstart"] = feimrch.nzstart;
602
603 dict["FEIMRC.mx"] = feimrch.mx;
604 dict["FEIMRC.my"] = feimrch.my;
605 dict["FEIMRC.mz"] = feimrch.mz;
606
607 dict["FEIMRC.xlen"] = feimrch.xlen;
608 dict["FEIMRC.ylen"] = feimrch.ylen;
609 dict["FEIMRC.zlen"] = feimrch.zlen;
610
611 dict["FEIMRC.alpha"] = feimrch.alpha;
612 dict["FEIMRC.beta"] = feimrch.beta;
613 dict["FEIMRC.gamma"] = feimrch.gamma;
614
615 dict["FEIMRC.mapc"] = feimrch.mapc;
616 dict["FEIMRC.mapr"] = feimrch.mapr;
617 dict["FEIMRC.maps"] = feimrch.maps;
618
619 dict["FEIMRC.minimum"] = feimrch.amin;
620 dict["FEIMRC.maximum"] = feimrch.amax;
621 dict["FEIMRC.mean"] = feimrch.amean;
622// dict["mean"] = feimrch.amean;
623
624 dict["FEIMRC.ispg"] = feimrch.ispg;
625 dict["FEIMRC.nsymbt"] = feimrch.nsymbt;
626
627 dict["apix_x"] = feimrch.xlen / feimrch.mx;
628 dict["apix_y"] = feimrch.ylen / feimrch.my;
629 dict["apix_z"] = feimrch.zlen / feimrch.mz;
630
631 dict["FEIMRC.next"] = feimrch.next; //offset from end of header to the first dataset
632 dict["FEIMRC.dvid"] = feimrch.dvid;
633 dict["FEIMRC.numintegers"] = feimrch.numintegers;
634 dict["FEIMRC.numfloats"] = feimrch.numfloats;
635 dict["FEIMRC.sub"] = feimrch.sub;
636 dict["FEIMRC.zfac"] = feimrch.zfac;
637
638 dict["FEIMRC.min2"] = feimrch.min2;
639 dict["FEIMRC.max2"] = feimrch.max2;
640 dict["FEIMRC.min3"] = feimrch.min3;
641 dict["FEIMRC.max3"] = feimrch.max3;
642 dict["FEIMRC.min4"] = feimrch.min4;
643 dict["FEIMRC.max4"] = feimrch.max4;
644
645 dict["FEIMRC.idtype"] = feimrch.idtype;
646 dict["FEIMRC.lens"] = feimrch.lens;
647 dict["FEIMRC.nd1"] = feimrch.nd1;
648 dict["FEIMRC.nd2"] = feimrch.nd2;
649 dict["FEIMRC.vd1"] = feimrch.vd1;
650 dict["FEIMRC.vd2"] = feimrch.vd2;
651
652 char label[32];
653
654 for(int i=0; i<9; i++) { // 9 tilt angles
655 sprintf(label, "MRC.tiltangles%d", i);
656 dict[string(label)] = feimrch.tiltangles[i];
657 }
658
659 dict["FEIMRC.zorg"] = feimrch.zorg;
660 dict["FEIMRC.xorg"] = feimrch.xorg;
661 dict["FEIMRC.yorg"] = feimrch.yorg;
662
663 dict["FEIMRC.nlabl"] = feimrch.nlabl;
664
665 for (int i = 0; i < feimrch.nlabl; i++) {
666 sprintf(label, "MRC.label%d", i);
667 dict[string(label)] = string(feimrch.labl[i], MRC_LABEL_SIZE);
668 }
669
670 /* Read extended image header by specified image index */
671 FeiMrcExtHeader feiexth;
672
673 portable_fseek(file, sizeof(FeiMrcHeader)+sizeof(FeiMrcExtHeader)*image_index, SEEK_SET);
674
675 if (fread(&feiexth, sizeof(FeiMrcExtHeader), 1, file) != 1) {
676 throw ImageReadException(filename, "FEI MRC extended header");
677 }
678
679 dict["FEIMRC.a_tilt"] = feiexth.a_tilt;
680 dict["FEIMRC.b_tilt"] = feiexth.b_tilt;
681
682 dict["FEIMRC.x_stage"] = feiexth.x_stage;
683 dict["FEIMRC.y_stage"] = feiexth.y_stage;
684 dict["FEIMRC.z_stage"] = feiexth.z_stage;
685
686 dict["FEIMRC.x_shift"] = feiexth.x_shift;
687 dict["FEIMRC.y_shift"] = feiexth.y_shift;
688
689 dict["FEIMRC.defocus"] = feiexth.defocus;
690 dict["FEIMRC.exp_time"] = feiexth.exp_time;
691 dict["FEIMRC.mean_int"] = feiexth.mean_int;
692 dict["FEIMRC.tilt_axis"] = feiexth.tilt_axis;
693
694 dict["FEIMRC.pixel_size"] = feiexth.pixel_size;
695 dict["apix_x"] = feiexth.pixel_size *1.0e10;
696 dict["apix_y"] = feiexth.pixel_size *1.0e10;
697 dict["apix_z"] = feiexth.pixel_size *1.0e10;
698 dict["FEIMRC.magnification"] = feiexth.magnification;
699 dict["FEIMRC.ht"] = feiexth.ht;
700 dict["FEIMRC.binning"] = feiexth.binning;
701 dict["FEIMRC.appliedDefocus"] = feiexth.appliedDefocus;
702
703 // remainder 16 4-byte floats not used
704
705 EXITFUNC;
706
707 return 0;
708}
709
710int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
711 EMUtil::EMDataType filestoragetype, bool use_host_endian)
712{
713 ENTERFUNC;
714
715 string ext = Util::get_filename_ext(filename);
716 is_stack = (ext == "mrcs" || ext == "MRCS");
717
718 bool append = (image_index == -1);
719
720 if (image_index == -1) {
721 image_index = 0;
722 }
723
724 if (image_index != 0 && ! is_stack) {
725 throw ImageWriteException(filename, "MRC file does not support stack.");
726 }
727
728 int max_images = 0;
729
730 if (! is_stack) {
731 max_images = 1;
732 }
733
734 check_write_access(rw_mode, image_index, max_images);
735
736 if (area) {
737 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
738
739 EXITFUNC;
740
741 return 0;
742 }
743
744 int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
745 int nx = dict["nx"];
746 int ny = dict["ny"];
747 int nz = dict["nz"];
748
749 bool got_one_image = (nz > 1);
750
751 is_ri = dict["is_complex_ri"];
752
753 bool opposite_endian = false;
754
755 if (! is_new_file) {
756 if (is_big_endian != ByteOrder::is_host_big_endian()) {
757 opposite_endian = true;
758 }
759
760 portable_fseek(file, 0, SEEK_SET);
761 }
762 else {
763 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
764 mrch.mapc = 1;
765 mrch.mapr = 2;
766 mrch.maps = 3;
767 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
768 }
769
770 if (nz <= 1 && dict.has_key("xform.projection") &&
771 ! dict.has_key("UCSF.chimera")) {
772 Transform * t = dict["xform.projection"];
773 Dict d = t->get_params("eman");
774// mrch.alpha = d["alpha"];
775// mrch.beta = d["beta"];
776// mrch.gamma = d["gamma"];
777 mrch.xorigin = d["tx"];
778 mrch.yorigin = d["ty"];
779 mrch.zorigin = d["tz"];
780
781 if (t) {delete t; t = NULL;}
782 }
783 else if (nz > 1 && dict.has_key("xform.align3d") &&
784 ! dict.has_key("UCSF.chimera")) {
785 Transform * t = dict["xform.align3d"];
786 Dict d = t->get_params("eman");
787// mrch.alpha = d["alpha"];
788// mrch.beta = d["beta"];
789// mrch.gamma = d["gamma"];
790 mrch.xorigin = d["tx"];
791 mrch.yorigin = d["ty"];
792 mrch.zorigin = d["tz"];
793
794 if (t) {delete t; t = NULL;}
795 }
796
797 if (dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
798 mrch.xorigin = (float)dict["origin_x"];
799 mrch.yorigin = (float)dict["origin_y"];
800
801 if (is_new_file) {
802 mrch.zorigin = (float)dict["origin_z"];
803 }
804 else {
805 mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
806 }
807 }
808
809 if (dict.has_key("MRC.nlabels")) {
810 mrch.nlabels = dict["MRC.nlabels"];
811 }
812
813 char label[32];
814
815 for (int i = 0; i < MRC_NUM_LABELS; i++) {
816 sprintf(label, "MRC.label%d", i);
817
818 if (dict.has_key(label)) {
819 strncpy(mrch.labels[i], (const char *) dict[label], MRC_LABEL_SIZE);
820 mrch.labels[i][MRC_LABEL_SIZE-1] = 0;
821
822 Util::replace_non_ascii(mrch.labels[i], MRC_LABEL_SIZE);
823
824 mrch.nlabels = i + 1;
825 }
826 }
827
828 if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
829 strcpy(mrch.labels[mrch.nlabels], "EMAN ");
830 strcat(mrch.labels[mrch.nlabels], Util::get_time_label().c_str());
831 mrch.nlabels++;
832 }
833
834 mrch.mode = new_mode;
835
836 if (is_complex_mode()) {
837 mrch.nx = nx / 2;
838 }
839 else {
840 mrch.nx = nx;
841 }
842
843 mrch.ny = ny;
844
845 if (is_stack) {
846 if (got_one_image) {
847 stack_size = nz;
848 image_index = 0;
849 }
850 else if (is_new_file) {
851 stack_size = 1;
852 image_index = stack_size - 1;
853 }
854 else if (append) {
855 stack_size++;
856 image_index = stack_size - 1;
857 }
858 else if (image_index >= stack_size) {
859 stack_size = image_index + 1;
860 }
861
862 nz = stack_size;
863 mrch.nz = nz;
864 }
865 else {
866 if (is_new_file) {
867 mrch.nz = nz;
868 }
869 else if (append) {
870 nz++;
871 mrch.nz = nz;
872 }
873 else if (image_index >= nz) {
874 nz = image_index + 1;
875 mrch.nz = nz;
876 }
877 else {
878 mrch.nz = nz;
879 }
880
881 stack_size = 1;
882 }
883
884 mrch.ispg = dict.has_key("MRC.ispg") ? (int)dict["MRC.ispg"] : 0;
885 mrch.nsymbt = 0;
886 mrch.amin = dict["minimum"];
887 mrch.amax = dict["maximum"];
888 mrch.amean = dict["mean"];
889 mrch.rms = dict["sigma"];
890
894// if(dict.has_key("MRC.mx")) {
895// mrch.mx = dict["MRC.mx"];
896// }
897// else {
898 mrch.mx = nx;
899// }
900// if(dict.has_key("MRC.my")) {
901// mrch.my = dict["MRC.my"];
902// }
903// else {
904 mrch.my = ny;
905// }
906// if(dict.has_key("MRC.mz")) {
907// mrch.mz = dict["MRC.mz"];
908// }
909// else {
910 mrch.mz = nz;
911// }
912
913 mrch.xlen = mrch.mx * (float) dict["apix_x"];
914 mrch.ylen = mrch.my * (float) dict["apix_y"];
915 mrch.zlen = mrch.mz * (float) dict["apix_z"];
916
917 if (dict.has_key("MRC.nxstart")) {
918 mrch.nxstart = dict["MRC.nxstart"];
919 }
920 else {
921 mrch.nxstart = -nx / 2;
922 }
923
924 if (dict.has_key("MRC.nystart")) {
925 mrch.nystart = dict["MRC.nystart"];
926 }
927 else {
928 mrch.nystart = -ny / 2;
929 }
930
931 if (dict.has_key("MRC.nzstart")) {
932 mrch.nzstart = dict["MRC.nzstart"];
933 }
934 else {
935 mrch.nzstart = -nz / 2;
936 }
937
938 strncpy(mrch.map, "MAP ", 4);
939 mrch.machinestamp = generate_machine_stamp();
940
941 MrcHeader mrch2 = mrch;
942
943 if (opposite_endian || !use_host_endian) {
944 swap_header(mrch2);
945 }
946
947 if (fwrite(&mrch2, sizeof(MrcHeader), 1, file) != 1) {
948 throw ImageWriteException(filename, "MRC header");
949 }
950
951 mode_size = get_mode_size(mrch.mode);
952 is_new_file = false;
953
954 if (is_stack && ! got_one_image) {
955 mrch.nz = 1;
956 }
957
958 EMUtil::getRenderLimits(dict, rendermin, rendermax, renderbits);
959
960 // Do not write ctf to mrc header in EMAN2
961
962// if( dict.has_key("ctf") ) {
963// vector<float> vctf = dict["ctf"];
964// EMAN1Ctf ctf_;
965// ctf_.from_vector(vctf);
966// write_ctf(ctf_);
967// }
968
969 EXITFUNC;
970
971 return 0;
972}
973
974int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool)
975{
976 ENTERFUNC;
977
978 if (! (isFEI || is_stack)) {
979 // single image format, index can only be zero
980 image_index = 0;
981 }
982
983 if (is_transpose && area != 0) {
984 printf("Warning: This image dimension is in (y,x,z), "
985 "region I/O not supported, return the whole image instead.");
986 }
987
988 check_read_access(image_index, rdata);
989
990 if (area && is_complex_mode()) {
991 LOGERR("Error: cannot read a region of a complex image.");
992
993 return 1;
994 }
995
996 signed char * scdata = (signed char *) rdata;
997 unsigned char * cdata = (unsigned char *) rdata;
998 short * sdata = (short *) rdata;
999 unsigned short * usdata = (unsigned short *) rdata;
1000
1001 size_t size = 0;
1002 int xlen = 0, ylen = 0, zlen = 0;
1003
1004 if (isFEI) { // FEI extended MRC
1005 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
1006 portable_fseek(file, sizeof(MrcHeader)+feimrch.next, SEEK_SET);
1007
1008 EMUtil::process_region_io(cdata, file, READ_ONLY,
1009 image_index, mode_size,
1010 feimrch.nx, feimrch.ny, feimrch.nz, area);
1011
1012 EMUtil::get_region_dims(area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
1013
1014 size = (size_t)xlen * ylen * zlen;
1015 }
1016 else { // regular MRC
1017 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
1018 portable_fseek(file, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
1019
1020 size_t modesize;
1021
1022 if (mrch.mode == MRC_UHEX) {
1023 // Have MRC packed 8 bit format with 2 4-bit values in each 8-bit byte,
1024 // so the mode size is effectively half a byte, signalled by this value:
1025 modesize = 11111111;
1026 }
1027 else {
1028 modesize = mode_size;
1029 }
1030
1031 EMUtil::process_region_io(cdata, file, READ_ONLY,
1032 image_index, modesize,
1033 mrch.nx, mrch.ny, mrch.nz, area);
1034
1035 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
1036
1037 size = (size_t)xlen * ylen * zlen;
1038 }
1039
1040 if (mrch.mode != MRC_UCHAR && mrch.mode != MRC_CHAR &&
1041 mrch.mode != MRC_UHEX) {
1042
1043 if (mode_size == sizeof(short)) {
1044 become_host_endian < short >(sdata, size);
1045 }
1046 else if (mode_size == sizeof(float)) {
1047 become_host_endian < float >(rdata, size);
1048 }
1049 }
1050
1051 if (mrch.mode == MRC_UHEX) {
1052 size_t num_pairs = size / 2;
1053 size_t num_pts = num_pairs * 2;
1054
1055 size_t ipt = num_pts;
1056
1057 for (size_t ipair = num_pairs; ipair >= 1; ipair--) {
1058 unsigned int v = (unsigned int) cdata[ipair-1];
1059 ipt--;
1060 rdata[ipt] = (float)(v >> 4); // v / 16;
1061 ipt--;
1062 rdata[ipt] = (float)(v & 15); // v % 16;
1063 }
1064 }
1065 else if (mrch.mode == MRC_UCHAR) {
1066 for (size_t i = 0; i < size; ++i) {
1067 size_t j = size - 1 - i;
1068 // rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
1069 rdata[j] = static_cast < float >(cdata[j]);
1070 }
1071 }
1072 else if (mrch.mode == MRC_CHAR) {
1073 for (size_t i = 0; i < size; ++i) {
1074 size_t j = size - 1 - i;
1075 // rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
1076 rdata[j] = static_cast < float >(scdata[j]);
1077 }
1078 }
1079 else if (mrch.mode == MRC_SHORT) {
1080 for (size_t i = 0; i < size; ++i) {
1081 size_t j = size - 1 - i;
1082 rdata[j] = static_cast < float >(sdata[j]);
1083 }
1084 }
1085 else if (mrch.mode == MRC_USHORT) {
1086 for (size_t i = 0; i < size; ++i) {
1087 size_t j = size - 1 - i;
1088 rdata[j] = static_cast < float >(usdata[j]);
1089 }
1090 }
1091
1092 if (is_transpose) {
1093 transpose(rdata, xlen, ylen, zlen);
1094 }
1095
1096 if (is_complex_mode()) {
1097 if (! is_ri) {
1098 Util::ap2ri(rdata, size);
1099 }
1100
1101 Util::flip_complex_phase(rdata, size);
1102 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
1103 }
1104
1105 EXITFUNC;
1106
1107 return 0;
1108}
1109
1110int MrcIO::write_data(float *data, int image_index, const Region* area,
1111 EMUtil::EMDataType, bool use_host_endian)
1112{
1113 ENTERFUNC;
1114
1115 bool append = (image_index == -1);
1116
1117 if (is_stack && append) {
1118 image_index = stack_size - 1;
1119 }
1120
1121 int max_images = 0;
1122
1123 if (! is_stack) {
1124 max_images = 1;
1125 image_index = 0;
1126 }
1127
1128 check_write_access(rw_mode, image_index, max_images, data);
1129 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
1130
1131 int nx, ny, nz;
1132
1133 if (! area) {
1134 nx = mrch.nx;
1135 ny = mrch.ny;
1136 nz = mrch.nz;
1137 }
1138 else {
1139 nx = (int)(area->get_width());
1140 ny = (int)(area->get_height());
1141 nz = (int)(area->get_depth());
1142 }
1143
1144 bool got_one_image = (nz > 1);
1145
1146 if (is_stack && ! got_one_image) {
1147 nz = 1;
1148 }
1149
1150 size_t size = (size_t)nx * ny * nz;
1151
1152 if (is_complex_mode()) {
1153 nx *= 2;
1154
1155 if (! is_ri) {
1156 Util::ap2ri(data, size);
1157 is_ri = 1;
1158 }
1159
1160 Util::flip_complex_phase(data, size);
1161 Util::rotate_phase_origin(data, nx, ny, nz);
1162 }
1163
1164 portable_fseek(file, sizeof(MrcHeader), SEEK_SET);
1165
1166 if ((is_big_endian != ByteOrder::is_host_big_endian()) || ! use_host_endian) {
1167 if (mrch.mode != MRC_UCHAR && mrch.mode != MRC_CHAR) {
1168 if (mode_size == sizeof(short)) {
1169 ByteOrder::swap_bytes((short*) data, size);
1170 }
1171 else if (mode_size == sizeof(float)) {
1172 ByteOrder::swap_bytes((float*) data, size);
1173 }
1174 }
1175 }
1176
1177 mode_size = get_mode_size(mrch.mode);
1178
1179// int xlen = 0, ylen = 0, zlen = 0;
1180// EMUtil::get_region_dims(area, nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
1181// int size = xlen * ylen * zlen;
1182 void * ptr_data = data;
1183
1184 int truebits=0;
1185 switch(mrch.mode) {
1186 case MRC_UCHAR: truebits=8; break;
1187 case MRC_SHORT: truebits=16; break;
1188 case MRC_FLOAT: truebits=32; break;
1189 case MRC_SHORT_COMPLEX: truebits=32; break;
1190 case MRC_FLOAT_COMPLEX: truebits=25; break;
1191 case MRC_USHORT: truebits=16; break;
1192 case MRC_UCHAR3: truebits=8; break;
1193 case MRC_CHAR: truebits=8; break;
1194 case MRC_UHEX: truebits=4; break;
1195 case MRC_UNKNOWN: truebits=0; break;
1196 }
1197 if (renderbits==0 || renderbits>truebits) renderbits=truebits;
1198
1199 float rmin = rendermin;
1200 float rmax = rendermax;
1201 int rbits = renderbits;
1202
1203 EMUtil::getRenderMinMax(data, nx, ny, rmin, rmax, rbits, nz);
1204
1205 signed char * scdata = NULL;
1206 unsigned char * cdata = NULL;
1207 short * sdata = NULL;
1208 unsigned short * usdata = NULL;
1209
1210 if (is_stack && ! got_one_image) {
1211 mrch.nz = stack_size;
1212 }
1213
1214 bool dont_rescale;
1215
1216 if (mrch.mode == MRC_UCHAR) {
1217 cdata = new unsigned char[size];
1218
1219 dont_rescale = (rmin == 0.0 && rmax == UCHAR_MAX);
1220
1221 for (size_t i = 0; i < size; ++i) {
1222 if (data[i] <= rmin) {
1223 cdata[i] = 0;
1224 }
1225 else if (data[i] >= rmax){
1226 cdata[i] = UCHAR_MAX;
1227 }
1228 else if (dont_rescale) {
1229 cdata[i] = (unsigned char) data[i];
1230 }
1231 else {
1232 cdata[i] = (unsigned char)((data[i] - rmin) /
1233 (rmax - rmin) *
1234 UCHAR_MAX);
1235 }
1236 }
1237
1238 ptr_data = cdata;
1239
1240 update_stats(ptr_data, size);
1241 }
1242 else if (mrch.mode == MRC_CHAR) {
1243 scdata = new signed char[size];
1244
1245 dont_rescale = (rmin == SCHAR_MIN && rmax == SCHAR_MAX);
1246
1247 for (size_t i = 0; i < size; ++i) {
1248 if (data[i] <= rmin) {
1249 scdata[i] = SCHAR_MIN;
1250 }
1251 else if (data[i] >= rmax){
1252 scdata[i] = SCHAR_MAX;
1253 }
1254 else if (dont_rescale) {
1255 scdata[i] = (signed char) data[i];
1256 }
1257 else {
1258 scdata[i] = (signed char)((data[i] - rmin) /
1259 (rmax - rmin) *
1260 (SCHAR_MAX - SCHAR_MIN) + SCHAR_MIN);
1261 }
1262 }
1263
1264 ptr_data = scdata;
1265
1266 update_stats(ptr_data, size);
1267 }
1268 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
1269 sdata = new short[size];
1270
1271 dont_rescale = (rmin == SHRT_MIN && rmax == SHRT_MAX);
1272
1273 for (size_t i = 0; i < size; ++i) {
1274 if (data[i] <= rmin) {
1275 sdata[i] = SHRT_MIN;
1276 }
1277 else if (data[i] >= rmax) {
1278 sdata[i] = SHRT_MAX;
1279 }
1280 else if (dont_rescale) {
1281 sdata[i] = (short) data[i];
1282 }
1283 else {
1284 sdata[i] = (short)(((data[i] - rmin) /
1285 (rmax - rmin)) *
1286 (SHRT_MAX - SHRT_MIN) + SHRT_MIN);
1287 }
1288 }
1289
1290 ptr_data = sdata;
1291
1292 update_stats(ptr_data, size);
1293 }
1294 else if (mrch.mode == MRC_USHORT) {
1295 usdata = new unsigned short[size];
1296
1297 dont_rescale = (rmin == 0.0 && rmax == USHRT_MAX);
1298
1299 for (size_t i = 0; i < size; ++i) {
1300 if (data[i] <= rmin) {
1301 usdata[i] = 0;
1302 }
1303 else if (data[i] >= rmax) {
1304 usdata[i] = USHRT_MAX;
1305 }
1306 else if (dont_rescale) {
1307 usdata[i] = (unsigned short) data[i];
1308 }
1309 else {
1310 usdata[i] = (unsigned short)((data[i] - rmin) /
1311 (rmax - rmin) *
1312 USHRT_MAX);
1313 }
1314 }
1315
1316 ptr_data = usdata;
1317
1318 update_stats(ptr_data, size);
1319 }
1320
1321 if (is_stack && ! got_one_image) {
1322 mrch.nz = 1;
1323 }
1324
1325 // New way to write data which includes region writing.
1326 // If it is tested to be OK, remove the old code in the
1327 // #if 0 ... #endif block.
1328 EMUtil::process_region_io(ptr_data, file, WRITE_ONLY, image_index,
1329 mode_size, mrch.nx, mrch.ny, mrch.nz, area);
1330
1331 if (cdata) {delete [] cdata; cdata = NULL;}
1332 if (scdata) {delete [] scdata; scdata = NULL;}
1333 if (sdata) {delete [] sdata; sdata = NULL;}
1334 if (usdata) {delete [] usdata; usdata = NULL;}
1335
1336#if 0
1337 int row_size = nx * get_mode_size(mrch.mode);
1338 int sec_size = nx * ny;
1339
1340 unsigned char * cbuf = new unsigned char[row_size];
1341 unsigned short * sbuf = (unsigned short *) cbuf;
1342
1343 for (int i = 0; i < nz; i++) {
1344 int i2 = i * sec_size;
1345
1346 for (int j = 0; j < ny; j++) {
1347 int k = i2 + j * nx;
1348 void * pbuf = 0;
1349
1350 switch (mrch.mode) {
1351 case MRC_UCHAR:
1352 for (int l = 0; l < nx; l++) {
1353 cbuf[l] = static_cast < unsigned char >(data[k + l]);
1354 }
1355
1356 pbuf = cbuf;
1357 fwrite(cbuf, row_size, 1, file);
1358
1359 break;
1360
1361 case MRC_SHORT:
1362 case MRC_SHORT_COMPLEX:
1363 for (int l = 0; l < nx; l++) {
1364 sbuf[l] = static_cast < short >(data[k + l]);
1365 }
1366
1367 pbuf = sbuf;
1368 fwrite(sbuf, row_size, 1, file);
1369
1370 break;
1371
1372 case MRC_USHORT:
1373 for (int l = 0; l < nx; l++) {
1374 sbuf[l] = static_cast < unsigned short >(data[k + l]);
1375 }
1376
1377 pbuf = sbuf;
1378 fwrite(sbuf, row_size, 1, file);
1379
1380 break;
1381
1382 case MRC_FLOAT:
1383 case MRC_FLOAT_COMPLEX:
1384 pbuf = &data[k];
1385
1386 break;
1387 }
1388
1389 if (pbuf) {
1390 fwrite(pbuf, row_size, 1, file);
1391 }
1392 }
1393 }
1394
1395 if (cbuf)
1396 {
1397 delete [] cbuf;
1398 cbuf = NULL;
1399 }
1400#endif
1401
1402 EXITFUNC;
1403 return 0;
1404}
1405
1406void MrcIO::update_stats(void * data, size_t size)
1407{
1408 float v; // variable to hold pixel value
1409 double sum;
1410 double square_sum;
1411 double mean;
1412 double sigma;
1413 double vv;
1414 float min, max;
1415
1416 signed char * scdata = NULL;
1417 unsigned char * cdata = NULL;
1418 short * sdata = NULL;
1419 unsigned short * usdata = NULL;
1420
1421 bool use_schar = (mrch.mode == MRC_CHAR);
1422 bool use_uchar = (mrch.mode == MRC_UCHAR);
1423 bool use_short = (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX);
1424 bool use_ushort = (mrch.mode == MRC_USHORT);
1425
1426 if (use_uchar) {
1427 max = 0.0;
1428 min = UCHAR_MAX;
1429 cdata = (unsigned char *) data;
1430 }
1431 else if (use_schar) {
1432 max = SCHAR_MIN;
1433 min = SCHAR_MAX;
1434 scdata = (signed char *) data;
1435 }
1436 else if (use_short) {
1437 max = (float) SHRT_MIN;
1438 min = (float) SHRT_MAX;
1439 sdata = (short *) data;
1440 }
1441 else if (use_ushort) {
1442 max = 0.0f;
1443 min = (float) USHRT_MAX;
1444 usdata = (unsigned short *) data;
1445 }
1446 else {
1447 throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
1448 }
1449
1450 sum = 0.0;
1451
1452 for (size_t i = 0; i < size; i++) {
1453 if (use_uchar) {
1454 v = (float) (cdata[i]);
1455 }
1456 else if (use_schar) {
1457 v = (float) (scdata[i]);
1458 }
1459 else if (use_short) {
1460 v = (float) (sdata[i]);
1461 }
1462 else {
1463 v = (float) (usdata[i]);
1464 }
1465
1466 if (v < min) min = v;
1467 if (v > max) max = v;
1468
1469 sum = sum + v;
1470 }
1471
1472 if (size > 0) {
1473 mean = sum / (double) size;
1474 }
1475 else {
1476 mean = 0.0;
1477 }
1478
1479 square_sum = 0.0;
1480
1481 for (size_t i = 0; i < size; i++) {
1482 if (use_uchar) {
1483 v = (float) (cdata[i]);
1484 }
1485 else if (use_schar) {
1486 v = (float) (scdata[i]);
1487 }
1488 else if (use_short) {
1489 v = (float) (sdata[i]);
1490 }
1491 else {
1492 v = (float) (usdata[i]);
1493 }
1494
1495 vv = v - mean;
1496
1497 square_sum = square_sum + vv * vv;
1498 }
1499
1500 if (size > 1) {
1501 sigma = std::sqrt(square_sum / (double) (size-1));
1502 }
1503 else {
1504 sigma = 0.0;
1505 }
1506
1507 /* change mrch.amin / amax / amean / rms here */
1508
1509 mrch.amin = min;
1510 mrch.amax = max;
1511 mrch.amean = (float) mean;
1512 mrch.rms = (float) sigma;
1513
1514// MrcHeader mrch2 = mrch;
1515//
1516// endian issue, can't get use_host_endian argument
1517// bool opposite_endian = false;
1518//
1519// if (!is_new_file) {
1520// if (is_big_endian != ByteOrder::is_host_big_endian()) {
1521// opposite_endian = true;
1522// }
1523//
1524// portable_fseek(mrcfile, 0, SEEK_SET);
1525// }
1526//
1527// if (opposite_endian || !use_host_endian) {
1528// swap_header(mrch2);
1529// }
1530
1531 portable_fseek(file, 0, SEEK_SET);
1532
1533 if (fwrite(& mrch, sizeof(MrcHeader), 1, file) != 1) {
1534 throw ImageWriteException(filename, "Error writing MRC header to update statistics.");
1535 }
1536
1537 portable_fseek(file, sizeof(MrcHeader), SEEK_SET);
1538}
1539
1540bool MrcIO::is_complex_mode()
1541{
1542 init();
1543
1544 if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
1545 return true;
1546 }
1547
1548 return false;
1549}
1550
1551int MrcIO::read_ctf(Ctf & ctf, int)
1552{
1553 ENTERFUNC;
1554
1555 init();
1556
1557 size_t n = strlen(CTF_MAGIC);
1558 int err = 1;
1559
1560 if (strncmp(mrch.labels[0], CTF_MAGIC, n) == 0) {
1561 err = ctf.from_string(string(&mrch.labels[0][n]));
1562 }
1563
1564 EXITFUNC;
1565
1566 return err;
1567}
1568
1569void MrcIO::write_ctf(const Ctf & ctf, int)
1570{
1571 ENTERFUNC;
1572
1573 init();
1574
1575 string ctf_str = ctf.to_string();
1576
1577 strcpy (mrch.labels[0], CTF_MAGIC);
1578 strncat(mrch.labels[0], ctf_str.c_str(),
1579 MRC_LABEL_SIZE - strlen(CTF_MAGIC) - 1);
1580
1581 rewind(file);
1582
1583 if (fwrite(&mrch, sizeof(MrcHeader), 1, file) != 1) {
1584 throw ImageWriteException(filename, "write CTF info to header failed");
1585 }
1586
1587 EXITFUNC;
1588}
1589
1590void MrcIO::flush()
1591{
1592 fflush(file);
1593}
1594
1595int MrcIO::get_mode_size(int mm)
1596{
1597 MrcIO::MrcMode m = static_cast < MrcMode > (mm);
1598
1599 int msize = 0;
1600 switch (m) {
1601 case MRC_CHAR:
1602 case MRC_UCHAR:
1603 msize = sizeof(char);
1604 break;
1605 case MRC_SHORT:
1606 case MRC_USHORT:
1607 case MRC_SHORT_COMPLEX:
1608 msize = sizeof(short);
1609 break;
1610 case MRC_FLOAT:
1611 case MRC_FLOAT_COMPLEX:
1612 msize = sizeof(float);
1613 break;
1614 default:
1615 msize = 0;
1616 }
1617
1618 return msize;
1619}
1620
1621int MrcIO::to_em_datatype(int m)
1622{
1623 EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
1624
1625 switch (m) {
1626 case MRC_CHAR:
1627 e = EMUtil::EM_CHAR;
1628 break;
1629 case MRC_UCHAR:
1630 e = EMUtil::EM_UCHAR;
1631 break;
1632 case MRC_SHORT:
1633 e = EMUtil::EM_SHORT;
1634 break;
1635 case MRC_USHORT:
1636 e = EMUtil::EM_USHORT;
1637 break;
1638 case MRC_SHORT_COMPLEX:
1639 e = EMUtil::EM_SHORT_COMPLEX;
1640 break;
1641 case MRC_FLOAT:
1642 e = EMUtil::EM_FLOAT;
1643 break;
1644 case MRC_FLOAT_COMPLEX:
1645 e = EMUtil::EM_FLOAT_COMPLEX;
1646 break;
1647 default:
1648 e = EMUtil::EM_UNKNOWN;
1649 }
1650
1651 return e;
1652}
1653
1654int MrcIO::to_mrcmode(int e, int is_complex)
1655{
1656 MrcMode m = MRC_UNKNOWN;
1657 EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
1658
1659 switch (em_type) {
1660 case EMUtil::EM_CHAR:
1661 m = MRC_CHAR;
1662
1663 break;
1664 case EMUtil::EM_UCHAR:
1665 m = MRC_UCHAR;
1666
1667 break;
1668 case EMUtil::EM_USHORT:
1669 if (is_complex) {
1670 m = MRC_SHORT_COMPLEX;
1671 }
1672 else {
1673 m = MRC_USHORT;
1674 }
1675
1676 break;
1677 case EMUtil::EM_SHORT:
1678 if (is_complex) {
1679 m = MRC_SHORT_COMPLEX;
1680 }
1681 else {
1682 m = MRC_SHORT;
1683 }
1684
1685 break;
1686 case EMUtil::EM_SHORT_COMPLEX:
1687 case EMUtil::EM_USHORT_COMPLEX:
1688 m = MRC_SHORT_COMPLEX;
1689
1690 break;
1691 case EMUtil::EM_INT:
1692 case EMUtil::EM_UINT:
1693 case EMUtil::EM_FLOAT:
1694 if (is_complex) {
1695 m = MRC_FLOAT_COMPLEX;
1696 }
1697 else {
1698 m = MRC_FLOAT;
1699 }
1700
1701 break;
1702 case EMUtil::EM_FLOAT_COMPLEX:
1703 m = MRC_FLOAT_COMPLEX;
1704
1705 break;
1706 default:
1707 m = MRC_FLOAT;
1708 }
1709
1710 return m;
1711}
1712
1713int MrcIO::generate_machine_stamp()
1714{
1715 int stamp = 0;
1716 char *p = (char *) (&stamp);
1717
1718 if (ByteOrder::is_host_big_endian()) {
1719 p[0] = 0x11;
1720 p[1] = 0x11;
1721 p[2] = 0;
1722 p[3] = 0;
1723 }
1724 else {
1725 p[0] = 0x44;
1726 p[1] = 0x41;
1727 p[2] = 0;
1728 p[3] = 0;
1729 }
1730
1731 return stamp;
1732}
1733
1734void MrcIO::swap_header(MrcHeader& mrch)
1735{
1736 ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
1737 ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
1738}
1739
1740int MrcIO::get_nimg()
1741{
1742 init();
1743
1744 return stack_size;
1745}
1746
1747int MrcIO::transpose(float *data, int xlen, int ylen, int zlen) const
1748{
1749 float * tmp = new float[xlen*ylen];
1750
1751 for(size_t z=0; z<(size_t)zlen; ++z) {
1752 for(size_t y=0; y<(size_t)ylen; ++y) {
1753 for(size_t x=0; x<(size_t)xlen; ++x) {
1754 tmp[x*ylen+y] = data[z*xlen*ylen+y*xlen+x];
1755 }
1756 }
1757
1758 std::copy(tmp, tmp+xlen*ylen, data+z*xlen*ylen);
1759 }
1760
1761 delete [] tmp;
1762
1763 return 0;
1764}
rdata
#define rdata(i)
Definition: analyzer.cpp:592
EMAN::ByteOrder::is_host_big_endian
static bool is_host_big_endian()
Definition: byteorder.cpp:40
EMAN::ByteOrder::swap_bytes
static void swap_bytes(T *data, size_t n=1)
swap the byte order of data with 'n' T-type elements.
Definition: byteorder.h:131
EMAN::Ctf
Ctf is the base class for all CTF model.
Definition: ctf.h:60
EMAN::Ctf::from_string
virtual int from_string(const string &ctf)=0
EMAN::Ctf::to_string
virtual string to_string() const =0
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::EMAN1Ctf
EMAN1Ctf is the CTF model used in EMAN1.
Definition: ctf.h:120
EMAN::EMAN1Ctf::to_vector
vector< float > to_vector() const
Definition: ctf.cpp:127
EMAN::EMUtil::process_region_io
static void process_region_io(void *cdata, FILE *file, int rw_mode, int image_index, size_t mode_size, int nx, int ny, int nz=1, const Region *area=0, bool need_flip=false, ImageType imgtype=IMAGE_UNKNOWN, int pre_row=0, int post_row=0)
Process image region IO.
Definition: emutil.cpp:931
EMAN::EMUtil::EMDataType
EMDataType
Image pixel data type used in EMAN.
Definition: emutil.h:92
EMAN::EMUtil::EM_UCHAR
@ EM_UCHAR
Definition: emutil.h:95
EMAN::EMUtil::EM_USHORT_COMPLEX
@ EM_USHORT_COMPLEX
Definition: emutil.h:103
EMAN::EMUtil::EM_UNKNOWN
@ EM_UNKNOWN
Definition: emutil.h:93
EMAN::EMUtil::EM_USHORT
@ EM_USHORT
Definition: emutil.h:97
EMAN::EMUtil::EM_FLOAT_COMPLEX
@ EM_FLOAT_COMPLEX
Definition: emutil.h:104
EMAN::EMUtil::EM_SHORT_COMPLEX
@ EM_SHORT_COMPLEX
Definition: emutil.h:102
EMAN::EMUtil::EM_SHORT
@ EM_SHORT
Definition: emutil.h:96
EMAN::EMUtil::getRenderMinMax
static void getRenderMinMax(float *data, const int nx, const int ny, float &rendermin, float &rendermax, int &renderbits, const int nz=1)
Calculate the min and max pixel value accepted for image nomalization, if we did not get them from im...
Definition: emutil.cpp:1706
EMAN::EMUtil::getRenderLimits
static void getRenderLimits(const Dict &dict, float &rendermin, float &rendermax, int &renderbits)
Get the min and max pixel value accepted for image nomalization from image attribute dictionary,...
Definition: emutil.cpp:1689
EMAN::EMUtil::get_region_dims
static void get_region_dims(const Region *area, int nx, int *area_x, int ny, int *area_y, int nz=1, int *area_z=0)
Get a region's dimensions.
Definition: emutil.cpp:860
EMAN::FloatSize
FloatSize is used to describe a 1D, 2D or 3D rectangular size in floating numbers.
Definition: geometry.h:105
EMAN::ImageIO
ImageIO classes are designed for reading/writing various electron micrography image formats,...
Definition: imageio.h:127
EMAN::ImageIO::rw_mode
IOMode rw_mode
Definition: imageio.h:353
EMAN::ImageIO::filename
string filename
Definition: imageio.h:352
EMAN::ImageIO::check_region
void check_region(const Region *area, const FloatSize &max_size, bool is_new_file=false, bool inbounds_only=true)
Validate image I/O region.
Definition: imageio.cpp:58
EMAN::ImageIO::flush
virtual void flush()=0
Flush the IO buffer.
EMAN::ImageIO::write_header
virtual int write_header(const Dict &dict, int image_index=0, const Region *area=0, EMUtil::EMDataType filestoragetype=EMUtil::EM_FLOAT, bool use_host_endian=true)=0
Write a header to an image.
EMAN::ImageIO::initialized
bool initialized
Definition: imageio.h:355
EMAN::ImageIO::write_data
virtual int write_data(float *data, int image_index=0, const Region *area=0, EMUtil::EMDataType filestoragetype=EMUtil::EM_FLOAT, bool use_host_endian=true)=0
Write data to an image.
EMAN::ImageIO::is_complex_mode
virtual bool is_complex_mode()=0
Is this an complex image or not.
EMAN::ImageIO::read_header
virtual int read_header(Dict &dict, int image_index=0, const Region *area=0, bool is_3d=false)=0
Read the header from an image.
EMAN::ImageIO::sfopen
FILE * sfopen(const string &filename, IOMode mode, bool *is_new=0, bool overwrite=false)
Run fopen safely.
Definition: imageio.cpp:135
EMAN::ImageIO::check_read_access
void check_read_access(int image_index)
Validate 'image_index' in file reading.
Definition: imageio.cpp:95
EMAN::ImageIO::read_data
virtual int read_data(float *data, int image_index=0, const Region *area=0, bool is_3d=false)=0
Read the data from an image.
EMAN::ImageIO::file
FILE * file
Definition: imageio.h:354
EMAN::ImageIO::init
virtual void init()=0
Do some initialization before doing the read/write.
EMAN::ImageIO::check_write_access
void check_write_access(IOMode rw_mode, int image_index, int max_nimg=0)
Validate rw_mode and image_index in file writing.
Definition: imageio.cpp:113
EMAN::ImageIO::is_image_big_endian
virtual bool is_image_big_endian()=0
Is this image in big endian or not.
EMAN::MrcIO::CTF_MAGIC
static const char * CTF_MAGIC
Definition: mrcio.h:260
EMAN::MrcIO::get_nimg
int get_nimg()
Return the number of images in this image file.
Definition: mrcio.cpp:1740
EMAN::MrcIO::check_swap
static void check_swap(const int *data, const char *filnam, bool show_errors, bool &do_swap, bool &have_err)
This is a utility routine to tell whether to byte swap MRC header.
Definition: mrcio.cpp:230
EMAN::MrcIO::read_fei_header
int read_fei_header(Dict &dict, int image_index=0, const Region *area=0, bool is_3d=false)
Definition: mrcio.cpp:574
EMAN::MrcIO::transpose
int transpose(float *data, int nx, int ny, int nz) const
Definition: mrcio.cpp:1747
EMAN::MrcIO::is_big_endian
bool is_big_endian
Definition: mrcio.h:283
EMAN::MrcIO::mode_size
int mode_size
Definition: mrcio.h:264
EMAN::MrcIO::rendermin
float rendermin
Definition: mrcio.h:286
EMAN::MrcIO::rendermax
float rendermax
Definition: mrcio.h:287
EMAN::MrcIO::is_valid
static bool is_valid(const void *first_block, off_t file_size=0)
Definition: mrcio.cpp:342
EMAN::MrcIO::renderbits
int renderbits
Definition: mrcio.h:288
EMAN::MrcIO::generate_machine_stamp
static int generate_machine_stamp()
generate the machine stamp used in MRC image format.
Definition: mrcio.cpp:1713
EMAN::MrcIO::use_given_dimensions
bool use_given_dimensions
Definition: mrcio.h:280
EMAN::MrcIO::write_ctf
void write_ctf(const Ctf &ctf, int image_index=0)
Write CTF data to this image.
Definition: mrcio.cpp:1569
EMAN::MrcIO::read_mrc_header
int read_mrc_header(Dict &dict, int image_index=0, const Region *area=0, bool is_3d=false)
Definition: mrcio.cpp:428
EMAN::MrcIO::to_em_datatype
static int to_em_datatype(int mrcmode)
Definition: mrcio.cpp:1621
EMAN::MrcIO::is_ri
int is_ri
Definition: mrcio.h:282
EMAN::MrcIO::MRC_LABEL_SIZE
@ MRC_LABEL_SIZE
Definition: mrcio.h:78
EMAN::MrcIO::NUM_4BYTES_AFTER_MAP
@ NUM_4BYTES_AFTER_MAP
Definition: mrcio.h:80
EMAN::MrcIO::NUM_4BYTES_PRE_MAP
@ NUM_4BYTES_PRE_MAP
Definition: mrcio.h:79
EMAN::MrcIO::MRC_NUM_LABELS
@ MRC_NUM_LABELS
Definition: mrcio.h:77
EMAN::MrcIO::is_stack
bool is_stack
Definition: mrcio.h:275
EMAN::MrcIO::is_transpose
bool is_transpose
Definition: mrcio.h:285
EMAN::MrcIO::mrch
MrcHeader mrch
Definition: mrcio.h:267
EMAN::MrcIO::swap_header
void swap_header(MrcHeader &mrch)
Definition: mrcio.cpp:1734
EMAN::MrcIO::isFEI
bool isFEI
Definition: mrcio.h:272
EMAN::MrcIO::get_mode_size
static int get_mode_size(int mm)
Definition: mrcio.cpp:1595
EMAN::MrcIO::is_8_bit_packed
bool is_8_bit_packed
Definition: mrcio.h:279
EMAN::MrcIO::feimrch
FeiMrcHeader feimrch
Definition: mrcio.h:268
EMAN::MrcIO::update_stats
void update_stats(void *data, size_t size)
This is a utility function used to calculate new min/max/mean/sigma when write MRC file as 16 bit or ...
Definition: mrcio.cpp:1406
EMAN::MrcIO::MRC_FLOAT_COMPLEX
@ MRC_FLOAT_COMPLEX
Definition: mrcio.h:68
EMAN::MrcIO::MRC_SHORT_COMPLEX
@ MRC_SHORT_COMPLEX
Definition: mrcio.h:67
EMAN::MrcIO::MRC_CHAR
@ MRC_CHAR
Definition: mrcio.h:71
EMAN::MrcIO::MRC_UCHAR3
@ MRC_UCHAR3
Definition: mrcio.h:70
EMAN::MrcIO::MRC_SHORT
@ MRC_SHORT
Definition: mrcio.h:65
EMAN::MrcIO::MRC_FLOAT
@ MRC_FLOAT
Definition: mrcio.h:66
EMAN::MrcIO::MRC_UNKNOWN
@ MRC_UNKNOWN
Definition: mrcio.h:73
EMAN::MrcIO::MRC_UCHAR
@ MRC_UCHAR
Definition: mrcio.h:64
EMAN::MrcIO::MRC_UHEX
@ MRC_UHEX
Definition: mrcio.h:72
EMAN::MrcIO::MRC_USHORT
@ MRC_USHORT
Definition: mrcio.h:69
EMAN::MrcIO::is_new_file
bool is_new_file
Definition: mrcio.h:284
EMAN::MrcIO::to_mrcmode
static int to_mrcmode(int em_datatype, int is_complex)
Definition: mrcio.cpp:1654
EMAN::MrcIO::stack_size
int stack_size
Definition: mrcio.h:276
EMAN::MrcIO::read_ctf
int read_ctf(Ctf &ctf, int image_index=0)
Read CTF data from this image.
Definition: mrcio.cpp:1551
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Region::get_width
float get_width() const
get the width
Definition: geometry.h:606
EMAN::Region::get_ndim
int get_ndim() const
Get the region's dimension.
Definition: geometry.h:644
EMAN::Region::get_depth
float get_depth() const
get the depth
Definition: geometry.h:610
EMAN::Region::origin
FloatPoint origin
Definition: geometry.h:654
EMAN::Region::get_height
float get_height() const
get the height
Definition: geometry.h:608
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::get_params
Dict get_params(const string &euler_type) const
Get the parameters of the entire transform, using a specific euler convention.
Definition: transform.cpp:479
EMAN::Util::flip_complex_phase
static void flip_complex_phase(float *data, size_t n)
flip the phase of a complex data array.
Definition: util.cpp:127
EMAN::Util::get_time_label
static string get_time_label()
Get the current time in a string with format "mm/dd/yyyy hh:mm".
Definition: util.cpp:1167
EMAN::Util::get_filename_ext
static string get_filename_ext(const string &filename)
Get a filename's extension.
Definition: util.cpp:526
EMAN::Util::replace_non_ascii
static void replace_non_ascii(char *str, int max_size, char repl_char='?')
Replace any non-ASCII characters in a C string with a given character.
Definition: util.cpp:289
EMAN::Util::rotate_phase_origin
static void rotate_phase_origin(float *data, size_t nx, size_t ny, size_t nz)
rotate data vertically by ny/2, to make the mrc phase origin compliance with EMAN1 and TNF reader
Definition: util.cpp:140
EMAN::Util::ap2ri
static void ap2ri(float *data, size_t n)
convert complex data array from Amplitude/Phase format into Real/Imaginary format.
Definition: util.cpp:97
EMAN::Util::get_min
static int get_min(int f1, int f2)
Get the minimum of 2 numbers.
Definition: util.h:922
copy
EMData * copy() const
This file is a part of "emdata.h", to use functions in this file, you should "#include "emdata....
sqrt
EMData * sqrt() const
return square root of current image
Definition: emdata_core.cpp:1084
is_ri
bool is_ri() const
Is this image a real/imaginary format complex image?
Definition: emdata_metadata.h:797
is_complex
bool is_complex() const
Is this a complex image?
Definition: emdata_metadata.h:648
ImageReadException
#define ImageReadException(filename, desc)
Definition: exception.h:204
ImageWriteException
#define ImageWriteException(imagename, desc)
Definition: exception.h:223
InvalidCallException
#define InvalidCallException(desc)
Definition: exception.h:348
LOGWARN
#define LOGWARN
Definition: log.h:53
LOGERR
#define LOGERR
Definition: log.h:51
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
EMAN
E2Exception class.
Definition: aligner.h:40
portable_fileio.h
portable_fseek
int portable_fseek(FILE *fp, off_t offset, int whence)
Definition: portable_fileio.h:39
y
#define y(i, j)
Definition: projector.cpp:1516
x
#define x(i)
Definition: projector.cpp:1517
EMAN::MrcIO::FeiMrcExtHeader
The extended header used by Fei MRC image.
Definition: mrcio.h:231
EMAN::MrcIO::FeiMrcHeader
Extended MRC format for tomography As used by Fei; original definition of extended header by Dave Aga...
Definition: mrcio.h:154
EMAN::MrcIO::FeiMrcHeader::labl
char labl[MRC_NUM_LABELS][MRC_LABEL_SIZE]
Definition: mrcio.h:222
EMAN::MrcIO::MrcHeader::labels
char labels[MRC_NUM_LABELS][MRC_LABEL_SIZE]
Definition: mrcio.h:138
Generated by doxygen 1.9.2