00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <cstring>
00037 #include <climits>
00038
00039 #include "mrcio.h"
00040 #include "portable_fileio.h"
00041 #include "geometry.h"
00042 #include "util.h"
00043 #include "ctf.h"
00044 #include "transform.h"
00045
00046 using namespace EMAN;
00047
00048 const char *MrcIO::CTF_MAGIC = "!-";
00049 const char *MrcIO::SHORT_CTF_MAGIC = "!$";
00050
00051 MrcIO::MrcIO(const string & mrc_filename, IOMode rw)
00052 : filename(mrc_filename), rw_mode(rw), mrcfile(0), mode_size(0)
00053 {
00054 memset(&mrch, 0, sizeof(MrcHeader));
00055 is_ri = 0;
00056 is_big_endian = ByteOrder::is_host_big_endian();
00057 is_new_file = false;
00058 initialized = false;
00059 }
00060
00061 MrcIO::~MrcIO()
00062 {
00063 if (mrcfile) {
00064 fclose(mrcfile);
00065 mrcfile = 0;
00066 }
00067 }
00068
00069 void MrcIO::init()
00070 {
00071 ENTERFUNC;
00072
00073 if (initialized) {
00074 return;
00075 }
00076
00077 initialized = true;
00078 mrcfile = sfopen(filename, rw_mode, &is_new_file);
00079
00080 if (!is_new_file) {
00081 if (fread(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00082 throw ImageReadException(filename, "MRC header");
00083 }
00084
00085 if (!is_valid(&mrch)) {
00086 throw ImageReadException(filename, "invalid MRC");
00087 }
00088
00089 is_big_endian = ByteOrder::is_data_big_endian(&mrch.nz);
00090 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00091 swap_header(mrch);
00092 }
00093
00094
00095 mode_size = get_mode_size(mrch.mode);
00096 if(is_complex_mode()) {
00097 is_ri = 1;
00098 }
00099
00100 if (mrch.nxstart != 0 || mrch.nystart != 0 || mrch.nzstart != 0) {
00101 LOGWARN("nx/ny/nz start not zero");
00102 }
00103
00104 if (is_complex_mode()) {
00105 mrch.nx *= 2;
00106 }
00107
00108 if (mrch.xlen == 0) {
00109 mrch.xlen = 1.0;
00110 }
00111
00112 if (mrch.ylen == 0) {
00113 mrch.ylen = 1.0;
00114 }
00115
00116 if (mrch.zlen == 0) {
00117 mrch.zlen = 1.0;
00118 }
00119 }
00120 EXITFUNC;
00121 }
00122
00123
00124 bool MrcIO::is_image_big_endian()
00125 {
00126 init();
00127 return is_big_endian;
00128 }
00129
00130 bool MrcIO::is_valid(const void *first_block, off_t file_size)
00131 {
00132 ENTERFUNC;
00133
00134 if (!first_block) {
00135 return false;
00136 }
00137
00138 const int *data = static_cast < const int *>(first_block);
00139 int nx = data[0];
00140 int ny = data[1];
00141 int nz = data[2];
00142 int mrcmode = data[3];
00143 int nsymbt = data[23];
00144
00145 bool data_big_endian = ByteOrder::is_data_big_endian(&nz);
00146
00147 if (data_big_endian != ByteOrder::is_host_big_endian()) {
00148 ByteOrder::swap_bytes(&nx);
00149 ByteOrder::swap_bytes(&ny);
00150 ByteOrder::swap_bytes(&nz);
00151 ByteOrder::swap_bytes(&mrcmode);
00152 }
00153
00154 if (mrcmode == MRC_SHORT_COMPLEX || mrcmode == MRC_FLOAT_COMPLEX) {
00155 nx *= 2;
00156 }
00157
00158 const int max_dim = 1 << 20;
00159
00160 if ((mrcmode >= MRC_UCHAR && mrcmode < MRC_UNKNOWN) &&
00161 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim) && (nz > 0 && nz < max_dim)) {
00162 #ifndef SPIDERMRC // Spider MRC files don't satisfy the following test
00163 if (file_size > 0) {
00164 off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz * (off_t)get_mode_size(mrcmode) + (off_t)sizeof(MrcHeader) + nsymbt;
00165 if (file_size == file_size1) {
00166 return true;
00167 }
00168 return false;
00169 }
00170 else {
00171 return true;
00172 }
00173 #endif // SPIDERMRC
00174 return true;
00175 }
00176 EXITFUNC;
00177 return false;
00178 }
00179
00180 int MrcIO::read_header(Dict & dict, int image_index, const Region * area, bool )
00181 {
00182 ENTERFUNC;
00183
00184
00185 if(image_index == -1) {
00186 image_index = 0;
00187 }
00188 if(image_index != 0) {
00189 throw ImageReadException(filename, "no stack allowed for MRC image. For take 2D slice out of 3D image, read the 3D image first, then use get_clip().");
00190 }
00191
00192 init();
00193
00194 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file,false);
00195
00196 dict["apix_x"] = mrch.xlen / mrch.mx;
00197 dict["apix_y"] = mrch.ylen / mrch.my;
00198 dict["apix_z"] = mrch.zlen / mrch.mz;
00199
00200 dict["MRC.minimum"] = mrch.amin;
00201 dict["MRC.maximum"] = mrch.amax;
00202 dict["MRC.mean"] = mrch.amean;
00203 dict["datatype"] = to_em_datatype(mrch.mode);
00204
00205 if (is_complex_mode()) {
00206 dict["is_complex"] = 1;
00207 dict["is_complex_ri"] = 1;
00208 }
00209
00210 int xlen = 0, ylen = 0, zlen = 0;
00211 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00212
00213 dict["nx"] = xlen;
00214 dict["ny"] = ylen;
00215 dict["nz"] = zlen;
00216
00217 if (area) {
00218 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
00219 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
00220
00221 if (area->get_ndim() == 3 && mrch.nz > 1) {
00222 dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
00223 }
00224 else {
00225 dict["origin_z"] = mrch.zorigin;
00226 }
00227 }
00228 else {
00229 dict["origin_x"] = mrch.xorigin;
00230 dict["origin_y"] = mrch.yorigin;
00231 dict["origin_z"] = mrch.zorigin;
00232 }
00233
00234 dict["MRC.nxstart"] = mrch.nxstart;
00235 dict["MRC.nystart"] = mrch.nystart;
00236 dict["MRC.nzstart"] = mrch.nzstart;
00237
00238 dict["MRC.mx"] = mrch.mx;
00239 dict["MRC.my"] = mrch.my;
00240 dict["MRC.mz"] = mrch.mz;
00241
00242 dict["MRC.nx"] = mrch.nx;
00243 dict["MRC.ny"] = mrch.ny;
00244 dict["MRC.nz"] = mrch.nz;
00245
00246 dict["MRC.xlen"] = mrch.xlen;
00247 dict["MRC.ylen"] = mrch.ylen;
00248 dict["MRC.zlen"] = mrch.zlen;
00249
00250 dict["MRC.alpha"] = mrch.alpha;
00251 dict["MRC.beta"] = mrch.beta;
00252 dict["MRC.gamma"] = mrch.gamma;
00253
00254 dict["MRC.mapc"] = mrch.mapc;
00255 dict["MRC.mapr"] = mrch.mapr;
00256 dict["MRC.maps"] = mrch.maps;
00257
00258 dict["MRC.ispg"] = mrch.ispg;
00259 dict["MRC.nsymbt"] = mrch.nsymbt;
00260 dict["MRC.machinestamp"] = mrch.machinestamp;
00261
00262 dict["MRC.rms"] = mrch.rms;
00263 dict["MRC.nlabels"] = mrch.nlabels;
00264 for (int i = 0; i < mrch.nlabels; i++) {
00265 char label[32];
00266 sprintf(label, "MRC.label%d", i);
00267 dict[string(label)] = mrch.labels[i];
00268 }
00269
00270 EMAN1Ctf ctf_;
00271 if(read_ctf(ctf_) == 0) {
00272 vector<float> vctf = ctf_.to_vector();
00273 dict["ctf"] = vctf;
00274 }
00275
00276 Dict dic;
00277 dic.put("type", "imagic");
00278 dic.put("alpha", mrch.alpha);
00279 dic.put("beta", mrch.beta);
00280 dic.put("gamma", mrch.gamma);
00281 dic.put("tx", mrch.xorigin);
00282 dic.put("ty", mrch.yorigin);
00283 dic.put("tz", mrch.zorigin);
00284 Transform * trans = new Transform(dic);
00285 if(zlen<=1) {
00286 dict["xform.projection"] = trans;
00287 }
00288 else {
00289 dict["xform.align3d"] = trans;
00290 }
00291
00292 if(trans) {delete trans; trans=0;}
00293 EXITFUNC;
00294 return 0;
00295 }
00296
00297 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00298 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00299 {
00300 ENTERFUNC;
00301
00302
00303 if(image_index == -1) {
00304 image_index = 0;
00305 }
00306 if(image_index != 0) {
00307 throw ImageWriteException(filename, "MRC file does not support stack.");
00308 }
00309 check_write_access(rw_mode, image_index, 1);
00310 if (area) {
00311 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00312 EXITFUNC;
00313 return 0;
00314 }
00315
00316 int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00317 int nx = dict["nx"];
00318 int ny = dict["ny"];
00319 int nz = dict["nz"];
00320 is_ri = dict["is_complex_ri"];
00321
00322 bool opposite_endian = false;
00323
00324 if (!is_new_file) {
00325 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00326 opposite_endian = true;
00327 }
00328 #if 0
00329 if (new_mode != mrch.mode) {
00330 LOGERR("cannot write to different mode file %s", filename.c_str());
00331 return 1;
00332 }
00333 #endif
00334 portable_fseek(mrcfile, 0, SEEK_SET);
00335 }
00336 else {
00337 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00338 mrch.mapc = 1;
00339 mrch.mapr = 2;
00340 mrch.maps = 3;
00341 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00342 }
00343
00344 if(nz<=1 && dict.has_key("xform.projection")) {
00345 Transform * t = dict["xform.projection"];
00346 Dict d = t->get_params("imagic");
00347 mrch.alpha = d["alpha"];
00348 mrch.beta = d["beta"];
00349 mrch.gamma = d["gamma"];
00350 mrch.xorigin = d["tx"];
00351 mrch.yorigin = d["ty"];
00352 mrch.zorigin = d["tz"];
00353 if(t) {delete t; t=0;}
00354 }
00355 else if(nz>1 && dict.has_key("xform.align3d")) {
00356 Transform * t = dict["xform.align3d"];
00357 Dict d = t->get_params("imagic");
00358 mrch.alpha = d["alpha"];
00359 mrch.beta = d["beta"];
00360 mrch.gamma = d["gamma"];
00361 mrch.xorigin = d["tx"];
00362 mrch.yorigin = d["ty"];
00363 mrch.zorigin = d["tz"];
00364 if(t) {delete t; t=0;}
00365 }
00366
00367 if(dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00368 mrch.xorigin = (float)dict["origin_x"];
00369 mrch.yorigin = (float)dict["origin_y"];
00370
00371 if (is_new_file) {
00372 mrch.zorigin = (float)dict["origin_z"];
00373 }
00374 else {
00375 mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00376 }
00377 }
00378
00379 if (dict.has_key("MRC.nlabels")) {
00380 mrch.nlabels = dict["MRC.nlabels"];
00381 }
00382
00383 for (int i = 0; i < MRC_NUM_LABELS; i++) {
00384 char label[32];
00385 sprintf(label, "MRC.label%d", i);
00386 if (dict.has_key(label)) {
00387 sprintf(&mrch.labels[i][0], "%s", (const char *) dict[label]);
00388 mrch.nlabels = i + 1;
00389 }
00390 }
00391
00392 if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00393 sprintf(&mrch.labels[mrch.nlabels][0], "EMAN %s", Util::get_time_label().c_str());
00394 mrch.nlabels++;
00395 }
00396
00397 mrch.labels[mrch.nlabels][0] = '\0';
00398 mrch.mode = new_mode;
00399
00400 if (is_complex_mode()) {
00401 mrch.nx = nx / 2;
00402 }
00403 else {
00404 mrch.nx = nx;
00405 }
00406 mrch.ny = ny;
00407
00408 if (is_new_file) {
00409 mrch.nz = nz;
00410 }
00411 else if (image_index >= mrch.nz) {
00412 mrch.nz = image_index + 1;
00413 }
00414
00415 mrch.ispg = 0;
00416 mrch.nsymbt = 0;
00417 mrch.amin = dict["minimum"];
00418 mrch.amax = dict["maximum"];
00419 mrch.amean = dict["mean"];
00420 mrch.rms = dict["sigma"];
00421
00424
00425
00426
00427
00428 mrch.mx = nx;
00429
00430
00431
00432
00433
00434 mrch.my = ny;
00435
00436
00437
00438
00439
00440 mrch.mz = nz;
00441
00442
00443 mrch.xlen = mrch.mx * (float) dict["apix_x"];
00444 mrch.ylen = mrch.my * (float) dict["apix_y"];
00445 mrch.zlen = mrch.mz * (float) dict["apix_z"];
00446
00447 if(dict.has_key("MRC.nxstart")) {
00448 mrch.nxstart = dict["MRC.nxstart"];
00449 }
00450 else {
00451 mrch.nxstart = -nx / 2;
00452 }
00453 if(dict.has_key("MRC.nystart")) {
00454 mrch.nystart = dict["MRC.nystart"];
00455 }
00456 else {
00457 mrch.nystart = -ny / 2;
00458 }
00459 if(dict.has_key("MRC.nzstart")) {
00460 mrch.nzstart = dict["MRC.nzstart"];
00461 }
00462 else {
00463 mrch.nzstart = -nz / 2;
00464 }
00465
00466 sprintf(mrch.map, "MAP ");
00467 mrch.machinestamp = generate_machine_stamp();
00468
00469 MrcHeader mrch2 = mrch;
00470
00471 if (opposite_endian || !use_host_endian) {
00472 swap_header(mrch2);
00473 }
00474
00475 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00476 throw ImageWriteException(filename, "MRC header");
00477 }
00478
00479 mode_size = get_mode_size(mrch.mode);
00480 is_new_file = false;
00481
00482 if( dict.has_key("ctf") ) {
00483 vector<float> vctf = dict["ctf"];
00484 EMAN1Ctf ctf_;
00485 ctf_.from_vector(vctf);
00486 write_ctf(ctf_);
00487 }
00488
00489 EXITFUNC;
00490 return 0;
00491 }
00492
00493 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool )
00494 {
00495 ENTERFUNC;
00496
00497
00498 image_index = 0;
00499 check_read_access(image_index, rdata);
00500
00501 if (area && is_complex_mode()) {
00502 LOGERR("Error: cannot read a region of a complex image.");
00503 return 1;
00504 }
00505
00506 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00507
00508 unsigned char *cdata = (unsigned char *) rdata;
00509 short *sdata = (short *) rdata;
00510 unsigned short *usdata = (unsigned short *) rdata;
00511
00512 portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00513
00514 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00515 image_index, mode_size,
00516 mrch.nx, mrch.ny, mrch.nz, area);
00517
00518 int xlen = 0, ylen = 0, zlen = 0;
00519 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00520
00521 size_t size = xlen * ylen * zlen;
00522
00523 if (mrch.mode != MRC_UCHAR) {
00524 if (mode_size == sizeof(short)) {
00525 become_host_endian < short >(sdata, size);
00526 }
00527 else if (mode_size == sizeof(float)) {
00528 become_host_endian < float >(rdata, size);
00529 }
00530 }
00531
00532 if (mrch.mode == MRC_UCHAR) {
00533 for (size_t i = 0; i < size; ++i) {
00534 size_t j = size - 1 - i;
00535
00536 rdata[j] = static_cast < float >(cdata[j]);
00537 }
00538 }
00539 else if (mrch.mode == MRC_SHORT ) {
00540 for (size_t i = 0; i < size; ++i) {
00541 size_t j = size - 1 - i;
00542 rdata[j] = static_cast < float >(sdata[j]);
00543 }
00544 }
00545 else if (mrch.mode == MRC_USHORT) {
00546 for (size_t i = 0; i < size; ++i) {
00547 size_t j = size - 1 - i;
00548 rdata[j] = static_cast < float >(usdata[j]);
00549 }
00550 }
00551
00552 if (is_complex_mode()) {
00553 if(!is_ri) Util::ap2ri(rdata, size);
00554 Util::flip_complex_phase(rdata, size);
00555 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00556 }
00557 EXITFUNC;
00558 return 0;
00559 }
00560
00561 int MrcIO::write_data(float *data, int image_index, const Region* area,
00562 EMUtil::EMDataType, bool use_host_endian)
00563 {
00564 ENTERFUNC;
00565
00566 image_index = 0;
00567 check_write_access(rw_mode, image_index, 1, data);
00568 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00569
00570 int nx = mrch.nx;
00571 int ny = mrch.ny;
00572 int nz = mrch.nz;
00573 size_t size = nx * ny * nz;
00574
00575 if (is_complex_mode()) {
00576 nx *= 2;
00577 if (!is_ri) {
00578 Util::ap2ri(data, size);
00579 is_ri = 1;
00580 }
00581 Util::flip_complex_phase(data, size);
00582 Util::rotate_phase_origin(data, nx, ny, nz);
00583 }
00584
00585 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00586
00587 if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian) {
00588 if (mrch.mode != MRC_UCHAR) {
00589 if (mode_size == sizeof(short)) {
00590 ByteOrder::swap_bytes((short*) data, size);
00591 }
00592 else if (mode_size == sizeof(float)) {
00593 ByteOrder::swap_bytes((float*) data, size);
00594 }
00595 }
00596 }
00597 mode_size = get_mode_size(mrch.mode);
00598
00599
00600
00601
00602 void * ptr_data = data;
00603
00604 float rendermin = 0.0f;
00605 float rendermax = 0.0f;
00606 getRenderMinMax(data, nx, ny, rendermin, rendermax);
00607
00608 unsigned char *cdata = 0;
00609 short *sdata = 0;
00610 unsigned short *usdata = 0;
00611 if (mrch.mode == MRC_UCHAR) {
00612 cdata = new unsigned char[size];
00613 for (size_t i = 0; i < size; ++i) {
00614 if(data[i] <= rendermin) {
00615 cdata[i] = 0;
00616 }
00617 else if(data[i] >= rendermax){
00618 cdata[i] = UCHAR_MAX;
00619 }
00620 else {
00621 cdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
00622 }
00623 }
00624 ptr_data = cdata;
00625 }
00626 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00627 sdata = new short[size];
00628 for (size_t i = 0; i < size; ++i) {
00629 if(data[i] <= rendermin) {
00630 sdata[i] = SHRT_MIN;
00631 }
00632 else if(data[i] >= rendermax) {
00633 sdata[i] = SHRT_MAX;
00634 }
00635 else {
00636 sdata[i]=(short)(((data[i]-rendermin)/(rendermax-rendermin))*(SHRT_MAX-SHRT_MIN) - SHRT_MAX);
00637 }
00638 }
00639 ptr_data = sdata;
00640 }
00641 else if (mrch.mode == MRC_USHORT) {
00642 usdata = new unsigned short[size];
00643 for (size_t i = 0; i < size; ++i) {
00644 if(data[i] <= rendermin) {
00645 usdata[i] = 0;
00646 }
00647 else if(data[i] >= rendermax) {
00648 usdata[i] = USHRT_MAX;
00649 }
00650 else {
00651 usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
00652 }
00653 }
00654 ptr_data = usdata;
00655 }
00656
00657
00658
00659
00660 EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
00661 mode_size, nx, mrch.ny, mrch.nz, area);
00662
00663 if(cdata) {delete [] cdata; cdata=0;}
00664 if(sdata) {delete [] sdata; sdata=0;}
00665 if(usdata) {delete [] usdata; usdata=0;}
00666
00667 #if 0
00668 int row_size = nx * get_mode_size(mrch.mode);
00669 int sec_size = nx * ny;
00670
00671 unsigned char *cbuf = new unsigned char[row_size];
00672 unsigned short *sbuf = (unsigned short *) cbuf;
00673
00674 for (int i = 0; i < nz; i++) {
00675 int i2 = i * sec_size;
00676 for (int j = 0; j < ny; j++) {
00677 int k = i2 + j * nx;
00678 void *pbuf = 0;
00679
00680 switch (mrch.mode) {
00681 case MRC_UCHAR:
00682 for (int l = 0; l < nx; l++) {
00683 cbuf[l] = static_cast < unsigned char >(data[k + l]);
00684 }
00685 pbuf = cbuf;
00686 fwrite(cbuf, row_size, 1, mrcfile);
00687 break;
00688
00689 case MRC_SHORT:
00690 case MRC_SHORT_COMPLEX:
00691 for (int l = 0; l < nx; l++) {
00692 sbuf[l] = static_cast < short >(data[k + l]);
00693 }
00694 pbuf = sbuf;
00695 fwrite(sbuf, row_size, 1, mrcfile);
00696 break;
00697
00698 case MRC_USHORT:
00699 for (int l = 0; l < nx; l++) {
00700 sbuf[l] = static_cast < unsigned short >(data[k + l]);
00701 }
00702 pbuf = sbuf;
00703 fwrite(sbuf, row_size, 1, mrcfile);
00704 break;
00705
00706 case MRC_FLOAT:
00707 case MRC_FLOAT_COMPLEX:
00708 pbuf = &data[k];
00709 break;
00710 }
00711 if (pbuf) {
00712 fwrite(pbuf, row_size, 1, mrcfile);
00713 }
00714 }
00715 }
00716
00717 if(cbuf)
00718 {
00719 delete[]cbuf;
00720 cbuf = 0;
00721 }
00722 #endif
00723
00724 EXITFUNC;
00725 return 0;
00726 }
00727
00728
00729 bool MrcIO::is_complex_mode()
00730 {
00731 init();
00732 if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
00733 return true;
00734 }
00735 return false;
00736 }
00737
00738
00739 int MrcIO::read_ctf(Ctf & ctf, int)
00740 {
00741 ENTERFUNC;
00742 init();
00743 size_t n = strlen(CTF_MAGIC);
00744
00745 int err = 1;
00746 if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
00747 err = ctf.from_string(string(&mrch.labels[0][n]));
00748 }
00749 EXITFUNC;
00750 return err;
00751 }
00752
00753 void MrcIO::write_ctf(const Ctf & ctf, int)
00754 {
00755 ENTERFUNC;
00756 init();
00757
00758 string ctf_str = ctf.to_string();
00759 sprintf(&mrch.labels[0][0], "%s%s", CTF_MAGIC, ctf_str.c_str());
00760 rewind(mrcfile);
00761
00762 if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00763 throw ImageWriteException(filename, "write CTF info to header failed");
00764 }
00765 EXITFUNC;
00766 }
00767
00768 void MrcIO::flush()
00769 {
00770 fflush(mrcfile);
00771 }
00772
00773
00774 int MrcIO::get_mode_size(int mm)
00775 {
00776 MrcIO::MrcMode m = static_cast < MrcMode > (mm);
00777
00778 int msize = 0;
00779 switch (m) {
00780 case MRC_UCHAR:
00781 msize = sizeof(char);
00782 break;
00783 case MRC_SHORT:
00784 case MRC_USHORT:
00785 case MRC_SHORT_COMPLEX:
00786 msize = sizeof(short);
00787 break;
00788 case MRC_FLOAT:
00789 case MRC_FLOAT_COMPLEX:
00790 msize = sizeof(float);
00791 break;
00792 default:
00793 msize = 0;
00794 }
00795
00796 return msize;
00797 }
00798
00799 int MrcIO::to_em_datatype(int m)
00800 {
00801 EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
00802
00803 switch (m) {
00804 case MRC_UCHAR:
00805 e = EMUtil::EM_UCHAR;
00806 break;
00807 case MRC_SHORT:
00808 e = EMUtil::EM_SHORT;
00809 break;
00810 case MRC_USHORT:
00811 e = EMUtil::EM_USHORT;
00812 break;
00813 case MRC_SHORT_COMPLEX:
00814 e = EMUtil::EM_SHORT_COMPLEX;
00815 break;
00816 case MRC_FLOAT:
00817 e = EMUtil::EM_FLOAT;
00818 break;
00819 case MRC_FLOAT_COMPLEX:
00820 e = EMUtil::EM_FLOAT_COMPLEX;
00821 break;
00822 default:
00823 e = EMUtil::EM_UNKNOWN;
00824 }
00825 return e;
00826 }
00827
00828
00829 int MrcIO::to_mrcmode(int e, int is_complex)
00830 {
00831 MrcMode m = MRC_UNKNOWN;
00832 EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
00833
00834 switch (em_type) {
00835 case EMUtil::EM_UCHAR:
00836 m = MRC_UCHAR;
00837 break;
00838 case EMUtil::EM_USHORT:
00839 if (is_complex) {
00840 m = MRC_SHORT_COMPLEX;
00841 }
00842 else {
00843 m = MRC_USHORT;
00844 }
00845 break;
00846 case EMUtil::EM_SHORT:
00847 if (is_complex) {
00848 m = MRC_SHORT_COMPLEX;
00849 }
00850 else {
00851 m = MRC_SHORT;
00852 }
00853 break;
00854 case EMUtil::EM_SHORT_COMPLEX:
00855 case EMUtil::EM_USHORT_COMPLEX:
00856 m = MRC_SHORT_COMPLEX;
00857 break;
00858 case EMUtil::EM_CHAR:
00859 case EMUtil::EM_INT:
00860 case EMUtil::EM_UINT:
00861 case EMUtil::EM_FLOAT:
00862 if (is_complex) {
00863 m = MRC_FLOAT_COMPLEX;
00864 }
00865 else {
00866 m = MRC_FLOAT;
00867 }
00868 break;
00869 case EMUtil::EM_FLOAT_COMPLEX:
00870 m = MRC_FLOAT_COMPLEX;
00871 break;
00872 default:
00873 m = MRC_FLOAT;
00874 }
00875
00876 return m;
00877 }
00878
00879
00880
00881 int MrcIO::generate_machine_stamp()
00882 {
00883 int stamp = 0;
00884 char *p = (char *) (&stamp);
00885
00886 if (ByteOrder::is_host_big_endian()) {
00887 p[0] = 0x11;
00888 p[1] = 0x11;
00889 p[2] = 0;
00890 p[3] = 0;
00891 }
00892 else {
00893 p[0] = 0x44;
00894 p[1] = 0x41;
00895 p[2] = 0;
00896 p[3] = 0;
00897 }
00898 return stamp;
00899 }
00900
00901 void MrcIO::swap_header(MrcHeader& mrch)
00902 {
00903 ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
00904 ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
00905 }