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 image_index = 0;
00186 check_read_access(image_index);
00187 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file,false);
00188
00189 dict["apix_x"] = mrch.xlen / mrch.mx;
00190 dict["apix_y"] = mrch.ylen / mrch.my;
00191 dict["apix_z"] = mrch.zlen / mrch.mz;
00192
00193 dict["minimum"] = mrch.amin;
00194 dict["maximum"] = mrch.amax;
00195 dict["mean"] = mrch.amean;
00196 dict["datatype"] = to_em_datatype(mrch.mode);
00197
00198 if (is_complex_mode()) {
00199 dict["is_complex"] = 1;
00200 dict["is_complex_ri"] = 1;
00201 }
00202
00203 int xlen = 0, ylen = 0, zlen = 0;
00204 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00205
00206 dict["nx"] = xlen;
00207 dict["ny"] = ylen;
00208 dict["nz"] = zlen;
00209
00210 if (area) {
00211 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
00212 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
00213
00214 if (area->get_ndim() == 3 && mrch.nz > 1) {
00215 dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
00216 }
00217 else {
00218 dict["origin_z"] = mrch.zorigin;
00219 }
00220 }
00221 else {
00222 dict["origin_x"] = mrch.xorigin;
00223 dict["origin_y"] = mrch.yorigin;
00224 dict["origin_z"] = mrch.zorigin;
00225 }
00226
00227 dict["MRC.nxstart"] = mrch.nxstart;
00228 dict["MRC.nystart"] = mrch.nystart;
00229 dict["MRC.nzstart"] = mrch.nzstart;
00230
00231 dict["MRC.mx"] = mrch.mx;
00232 dict["MRC.my"] = mrch.my;
00233 dict["MRC.mz"] = mrch.mz;
00234
00235 dict["MRC.nx"] = mrch.nx;
00236 dict["MRC.ny"] = mrch.ny;
00237 dict["MRC.nz"] = mrch.nz;
00238
00239 dict["MRC.xlen"] = mrch.xlen;
00240 dict["MRC.ylen"] = mrch.ylen;
00241 dict["MRC.zlen"] = mrch.zlen;
00242
00243 dict["MRC.alpha"] = mrch.alpha;
00244 dict["MRC.beta"] = mrch.beta;
00245 dict["MRC.gamma"] = mrch.gamma;
00246
00247 dict["MRC.mapc"] = mrch.mapc;
00248 dict["MRC.mapr"] = mrch.mapr;
00249 dict["MRC.maps"] = mrch.maps;
00250
00251 dict["MRC.ispg"] = mrch.ispg;
00252 dict["MRC.nsymbt"] = mrch.nsymbt;
00253 dict["MRC.machinestamp"] = mrch.machinestamp;
00254
00255 dict["MRC.rms"] = mrch.rms;
00256 dict["MRC.nlabels"] = mrch.nlabels;
00257 for (int i = 0; i < mrch.nlabels; i++) {
00258 char label[32];
00259 sprintf(label, "MRC.label%d", i);
00260 dict[string(label)] = mrch.labels[i];
00261 }
00262
00263 EMAN1Ctf ctf_;
00264 if(read_ctf(ctf_) == 0) {
00265 vector<float> vctf = ctf_.to_vector();
00266 dict["ctf"] = vctf;
00267 }
00268
00269 Dict dic;
00270 dic.put("type", "imagic");
00271 dic.put("alpha", mrch.alpha);
00272 dic.put("beta", mrch.beta);
00273 dic.put("gamma", mrch.gamma);
00274 dic.put("tx", mrch.xorigin);
00275 dic.put("ty", mrch.yorigin);
00276 dic.put("tz", mrch.zorigin);
00277 Transform * trans = new Transform(dic);
00278 if(zlen<=1) {
00279 dict["xform.projection"] = trans;
00280 }
00281 else {
00282 dict["xform.align3d"] = trans;
00283 }
00284
00285 if(trans) {delete trans; trans=0;}
00286 EXITFUNC;
00287 return 0;
00288 }
00289
00290 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00291 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00292 {
00293 ENTERFUNC;
00294
00295
00296 image_index = 0;
00297 check_write_access(rw_mode, image_index, 1);
00298 if (area) {
00299 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00300 EXITFUNC;
00301 return 0;
00302 }
00303
00304 int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00305 int nx = dict["nx"];
00306 int ny = dict["ny"];
00307 int nz = dict["nz"];
00308 is_ri = dict["is_complex_ri"];
00309
00310 bool opposite_endian = false;
00311
00312 if (!is_new_file) {
00313 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00314 opposite_endian = true;
00315 }
00316 #if 0
00317 if (new_mode != mrch.mode) {
00318 LOGERR("cannot write to different mode file %s", filename.c_str());
00319 return 1;
00320 }
00321 #endif
00322 portable_fseek(mrcfile, 0, SEEK_SET);
00323 }
00324 else {
00325 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00326 mrch.mapc = 1;
00327 mrch.mapr = 2;
00328 mrch.maps = 3;
00329 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00330 }
00331
00332 if(nz<=1 && dict.has_key("xform.projection")) {
00333 Transform * t = dict["xform.projection"];
00334 Dict d = t->get_params("imagic");
00335 mrch.alpha = d["alpha"];
00336 mrch.beta = d["beta"];
00337 mrch.gamma = d["gamma"];
00338 mrch.xorigin = d["tx"];
00339 mrch.yorigin = d["ty"];
00340 mrch.zorigin = d["tz"];
00341 if(t) {delete t; t=0;}
00342 }
00343 else if(nz>1 && dict.has_key("xform.align3d")) {
00344 Transform * t = dict["xform.align3d"];
00345 Dict d = t->get_params("imagic");
00346 mrch.alpha = d["alpha"];
00347 mrch.beta = d["beta"];
00348 mrch.gamma = d["gamma"];
00349 mrch.xorigin = d["tx"];
00350 mrch.yorigin = d["ty"];
00351 mrch.zorigin = d["tz"];
00352 if(t) {delete t; t=0;}
00353 }
00354
00355 if(dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00356 mrch.xorigin = (float)dict["origin_x"];
00357 mrch.yorigin = (float)dict["origin_y"];
00358
00359 if (is_new_file) {
00360 mrch.zorigin = (float)dict["origin_z"];
00361 }
00362 else {
00363 mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00364 }
00365 }
00366
00367 if (dict.has_key("MRC.nlabels")) {
00368 mrch.nlabels = dict["MRC.nlabels"];
00369 }
00370
00371 for (int i = 0; i < MRC_NUM_LABELS; i++) {
00372 char label[32];
00373 sprintf(label, "MRC.label%d", i);
00374 if (dict.has_key(label)) {
00375 sprintf(&mrch.labels[i][0], "%s", (const char *) dict[label]);
00376 mrch.nlabels = i + 1;
00377 }
00378 }
00379
00380 if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00381 sprintf(&mrch.labels[mrch.nlabels][0], "EMAN %s", Util::get_time_label().c_str());
00382 mrch.nlabels++;
00383 }
00384
00385 mrch.labels[mrch.nlabels][0] = '\0';
00386 mrch.mode = new_mode;
00387
00388 if (is_complex_mode()) {
00389 mrch.nx = nx / 2;
00390 }
00391 else {
00392 mrch.nx = nx;
00393 }
00394 mrch.ny = ny;
00395
00396 if (is_new_file) {
00397 mrch.nz = nz;
00398 }
00399 else if (image_index >= mrch.nz) {
00400 mrch.nz = image_index + 1;
00401 }
00402
00403 mrch.ispg = 0;
00404 mrch.nsymbt = 0;
00405 mrch.amin = dict["minimum"];
00406 mrch.amax = dict["maximum"];
00407 mrch.amean = dict["mean"];
00408
00411
00412
00413
00414
00415 mrch.mx = nx;
00416
00417
00418
00419
00420
00421 mrch.my = ny;
00422
00423
00424
00425
00426
00427 mrch.mz = nz;
00428
00429
00430 mrch.xlen = mrch.mx * (float) dict["apix_x"];
00431 mrch.ylen = mrch.my * (float) dict["apix_y"];
00432 mrch.zlen = mrch.mz * (float) dict["apix_z"];
00433
00434 if(dict.has_key("MRC.nxstart")) {
00435 mrch.nxstart = dict["MRC.nxstart"];
00436 }
00437 else {
00438 mrch.nxstart = -nx / 2;
00439 }
00440 if(dict.has_key("MRC.nystart")) {
00441 mrch.nystart = dict["MRC.nystart"];
00442 }
00443 else {
00444 mrch.nystart = -ny / 2;
00445 }
00446 if(dict.has_key("MRC.nzstart")) {
00447 mrch.nzstart = dict["MRC.nzstart"];
00448 }
00449 else {
00450 mrch.nzstart = -nz / 2;
00451 }
00452
00453 sprintf(mrch.map, "MAP ");
00454 mrch.machinestamp = generate_machine_stamp();
00455 if(dict.has_key("MRC.rms")) {
00456 mrch.rms = (float)dict["MRC.rms"];
00457 }
00458
00459 MrcHeader mrch2 = mrch;
00460
00461 if (opposite_endian || !use_host_endian) {
00462 swap_header(mrch2);
00463 }
00464
00465 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00466 throw ImageWriteException(filename, "MRC header");
00467 }
00468
00469 mode_size = get_mode_size(mrch.mode);
00470 is_new_file = false;
00471
00472 if( dict.has_key("ctf") ) {
00473 vector<float> vctf = dict["ctf"];
00474 EMAN1Ctf ctf_;
00475 ctf_.from_vector(vctf);
00476 write_ctf(ctf_);
00477 }
00478
00479 EXITFUNC;
00480 return 0;
00481 }
00482
00483 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool )
00484 {
00485 ENTERFUNC;
00486
00487
00488 image_index = 0;
00489 check_read_access(image_index, rdata);
00490
00491 if (area && is_complex_mode()) {
00492 LOGERR("Error: cannot read a region of a complex image.");
00493 return 1;
00494 }
00495
00496 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00497
00498 unsigned char *cdata = (unsigned char *) rdata;
00499 short *sdata = (short *) rdata;
00500 unsigned short *usdata = (unsigned short *) rdata;
00501
00502 portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00503
00504 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00505 image_index, mode_size,
00506 mrch.nx, mrch.ny, mrch.nz, area);
00507
00508 int xlen = 0, ylen = 0, zlen = 0;
00509 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00510
00511 size_t size = xlen * ylen * zlen;
00512
00513 if (mrch.mode != MRC_UCHAR) {
00514 if (mode_size == sizeof(short)) {
00515 become_host_endian < short >(sdata, size);
00516 }
00517 else if (mode_size == sizeof(float)) {
00518 become_host_endian < float >(rdata, size);
00519 }
00520 }
00521
00522 if (mrch.mode == MRC_UCHAR) {
00523 for (size_t i = 0; i < size; ++i) {
00524 size_t j = size - 1 - i;
00525
00526 rdata[j] = static_cast < float >(cdata[j]);
00527 }
00528 }
00529 else if (mrch.mode == MRC_SHORT ) {
00530 for (size_t i = 0; i < size; ++i) {
00531 size_t j = size - 1 - i;
00532 rdata[j] = static_cast < float >(sdata[j]);
00533 }
00534 }
00535 else if (mrch.mode == MRC_USHORT) {
00536 for (size_t i = 0; i < size; ++i) {
00537 size_t j = size - 1 - i;
00538 rdata[j] = static_cast < float >(usdata[j]);
00539 }
00540 }
00541
00542 if (is_complex_mode()) {
00543 if(!is_ri) Util::ap2ri(rdata, size);
00544 Util::flip_complex_phase(rdata, size);
00545 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00546 }
00547 EXITFUNC;
00548 return 0;
00549 }
00550
00551 int MrcIO::write_data(float *data, int image_index, const Region* area,
00552 EMUtil::EMDataType, bool use_host_endian)
00553 {
00554 ENTERFUNC;
00555
00556 image_index = 0;
00557 check_write_access(rw_mode, image_index, 1, data);
00558 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00559
00560 int nx = mrch.nx;
00561 int ny = mrch.ny;
00562 int nz = mrch.nz;
00563 size_t size = nx * ny * nz;
00564
00565 if (is_complex_mode()) {
00566 nx *= 2;
00567 if (!is_ri) {
00568 Util::ap2ri(data, size);
00569 is_ri = 1;
00570 }
00571 Util::flip_complex_phase(data, size);
00572 Util::rotate_phase_origin(data, nx, ny, nz);
00573 }
00574
00575 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00576
00577 if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian) {
00578 if (mrch.mode != MRC_UCHAR) {
00579 if (mode_size == sizeof(short)) {
00580 ByteOrder::swap_bytes((short*) data, size);
00581 }
00582 else if (mode_size == sizeof(float)) {
00583 ByteOrder::swap_bytes((float*) data, size);
00584 }
00585 }
00586 }
00587 mode_size = get_mode_size(mrch.mode);
00588
00589
00590
00591
00592 void * ptr_data = data;
00593
00594 float rendermin = 0.0f;
00595 float rendermax = 0.0f;
00596 getRenderMinMax(data, nx, ny, rendermin, rendermax);
00597
00598 unsigned char *cdata = 0;
00599 short *sdata = 0;
00600 unsigned short *usdata = 0;
00601 if (mrch.mode == MRC_UCHAR) {
00602 cdata = new unsigned char[size];
00603 for (size_t i = 0; i < size; ++i) {
00604 if(data[i] <= rendermin) {
00605 cdata[i] = 0;
00606 }
00607 else if(data[i] >= rendermax){
00608 cdata[i] = UCHAR_MAX;
00609 }
00610 else {
00611 cdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
00612 }
00613 }
00614 ptr_data = cdata;
00615 }
00616 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00617 sdata = new short[size];
00618 for (size_t i = 0; i < size; ++i) {
00619 if(data[i] <= rendermin) {
00620 sdata[i] = SHRT_MIN;
00621 }
00622 else if(data[i] >= rendermax) {
00623 sdata[i] = SHRT_MAX;
00624 }
00625 else {
00626 sdata[i]=(short)(((data[i]-rendermin)/(rendermax-rendermin))*(SHRT_MAX-SHRT_MIN) - SHRT_MAX);
00627 }
00628 }
00629 ptr_data = sdata;
00630 }
00631 else if (mrch.mode == MRC_USHORT) {
00632 usdata = new unsigned short[size];
00633 for (size_t i = 0; i < size; ++i) {
00634 if(data[i] <= rendermin) {
00635 usdata[i] = 0;
00636 }
00637 else if(data[i] >= rendermax) {
00638 usdata[i] = USHRT_MAX;
00639 }
00640 else {
00641 usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
00642 }
00643 }
00644 ptr_data = usdata;
00645 }
00646
00647
00648
00649
00650 EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
00651 mode_size, nx, mrch.ny, mrch.nz, area);
00652
00653 if(cdata) {delete [] cdata; cdata=0;}
00654 if(sdata) {delete [] sdata; sdata=0;}
00655 if(usdata) {delete [] usdata; usdata=0;}
00656
00657 #if 0
00658 int row_size = nx * get_mode_size(mrch.mode);
00659 int sec_size = nx * ny;
00660
00661 unsigned char *cbuf = new unsigned char[row_size];
00662 unsigned short *sbuf = (unsigned short *) cbuf;
00663
00664 for (int i = 0; i < nz; i++) {
00665 int i2 = i * sec_size;
00666 for (int j = 0; j < ny; j++) {
00667 int k = i2 + j * nx;
00668 void *pbuf = 0;
00669
00670 switch (mrch.mode) {
00671 case MRC_UCHAR:
00672 for (int l = 0; l < nx; l++) {
00673 cbuf[l] = static_cast < unsigned char >(data[k + l]);
00674 }
00675 pbuf = cbuf;
00676 fwrite(cbuf, row_size, 1, mrcfile);
00677 break;
00678
00679 case MRC_SHORT:
00680 case MRC_SHORT_COMPLEX:
00681 for (int l = 0; l < nx; l++) {
00682 sbuf[l] = static_cast < short >(data[k + l]);
00683 }
00684 pbuf = sbuf;
00685 fwrite(sbuf, row_size, 1, mrcfile);
00686 break;
00687
00688 case MRC_USHORT:
00689 for (int l = 0; l < nx; l++) {
00690 sbuf[l] = static_cast < unsigned short >(data[k + l]);
00691 }
00692 pbuf = sbuf;
00693 fwrite(sbuf, row_size, 1, mrcfile);
00694 break;
00695
00696 case MRC_FLOAT:
00697 case MRC_FLOAT_COMPLEX:
00698 pbuf = &data[k];
00699 break;
00700 }
00701 if (pbuf) {
00702 fwrite(pbuf, row_size, 1, mrcfile);
00703 }
00704 }
00705 }
00706
00707 if(cbuf)
00708 {
00709 delete[]cbuf;
00710 cbuf = 0;
00711 }
00712 #endif
00713
00714 EXITFUNC;
00715 return 0;
00716 }
00717
00718
00719 bool MrcIO::is_complex_mode()
00720 {
00721 init();
00722 if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
00723 return true;
00724 }
00725 return false;
00726 }
00727
00728
00729 int MrcIO::read_ctf(Ctf & ctf, int)
00730 {
00731 ENTERFUNC;
00732 init();
00733 size_t n = strlen(CTF_MAGIC);
00734
00735 int err = 1;
00736 if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
00737 err = ctf.from_string(string(&mrch.labels[0][n]));
00738 }
00739 EXITFUNC;
00740 return err;
00741 }
00742
00743 void MrcIO::write_ctf(const Ctf & ctf, int)
00744 {
00745 ENTERFUNC;
00746 init();
00747
00748 string ctf_str = ctf.to_string();
00749 sprintf(&mrch.labels[0][0], "%s%s", CTF_MAGIC, ctf_str.c_str());
00750 rewind(mrcfile);
00751
00752 if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00753 throw ImageWriteException(filename, "write CTF info to header failed");
00754 }
00755 EXITFUNC;
00756 }
00757
00758 void MrcIO::flush()
00759 {
00760 fflush(mrcfile);
00761 }
00762
00763
00764 int MrcIO::get_mode_size(int mm)
00765 {
00766 MrcIO::MrcMode m = static_cast < MrcMode > (mm);
00767
00768 int msize = 0;
00769 switch (m) {
00770 case MRC_UCHAR:
00771 msize = sizeof(char);
00772 break;
00773 case MRC_SHORT:
00774 case MRC_USHORT:
00775 case MRC_SHORT_COMPLEX:
00776 msize = sizeof(short);
00777 break;
00778 case MRC_FLOAT:
00779 case MRC_FLOAT_COMPLEX:
00780 msize = sizeof(float);
00781 break;
00782 default:
00783 msize = 0;
00784 }
00785
00786 return msize;
00787 }
00788
00789 int MrcIO::to_em_datatype(int m)
00790 {
00791 EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
00792
00793 switch (m) {
00794 case MRC_UCHAR:
00795 e = EMUtil::EM_UCHAR;
00796 break;
00797 case MRC_SHORT:
00798 e = EMUtil::EM_SHORT;
00799 break;
00800 case MRC_USHORT:
00801 e = EMUtil::EM_USHORT;
00802 break;
00803 case MRC_SHORT_COMPLEX:
00804 e = EMUtil::EM_SHORT_COMPLEX;
00805 break;
00806 case MRC_FLOAT:
00807 e = EMUtil::EM_FLOAT;
00808 break;
00809 case MRC_FLOAT_COMPLEX:
00810 e = EMUtil::EM_FLOAT_COMPLEX;
00811 break;
00812 default:
00813 e = EMUtil::EM_UNKNOWN;
00814 }
00815 return e;
00816 }
00817
00818
00819 int MrcIO::to_mrcmode(int e, int is_complex)
00820 {
00821 MrcMode m = MRC_UNKNOWN;
00822 EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
00823
00824 switch (em_type) {
00825 case EMUtil::EM_UCHAR:
00826 m = MRC_UCHAR;
00827 break;
00828 case EMUtil::EM_USHORT:
00829 if (is_complex) {
00830 m = MRC_SHORT_COMPLEX;
00831 }
00832 else {
00833 m = MRC_USHORT;
00834 }
00835 break;
00836 case EMUtil::EM_SHORT:
00837 if (is_complex) {
00838 m = MRC_SHORT_COMPLEX;
00839 }
00840 else {
00841 m = MRC_SHORT;
00842 }
00843 break;
00844 case EMUtil::EM_SHORT_COMPLEX:
00845 case EMUtil::EM_USHORT_COMPLEX:
00846 m = MRC_SHORT_COMPLEX;
00847 break;
00848 case EMUtil::EM_CHAR:
00849 case EMUtil::EM_INT:
00850 case EMUtil::EM_UINT:
00851 case EMUtil::EM_FLOAT:
00852 if (is_complex) {
00853 m = MRC_FLOAT_COMPLEX;
00854 }
00855 else {
00856 m = MRC_FLOAT;
00857 }
00858 break;
00859 case EMUtil::EM_FLOAT_COMPLEX:
00860 m = MRC_FLOAT_COMPLEX;
00861 break;
00862 default:
00863 m = MRC_FLOAT;
00864 }
00865
00866 return m;
00867 }
00868
00869
00870
00871 int MrcIO::generate_machine_stamp()
00872 {
00873 int stamp = 0;
00874 char *p = (char *) (&stamp);
00875
00876 if (ByteOrder::is_host_big_endian()) {
00877 p[0] = 0x11;
00878 p[1] = 0x11;
00879 p[2] = 0;
00880 p[3] = 0;
00881 }
00882 else {
00883 p[0] = 0x44;
00884 p[1] = 0x41;
00885 p[2] = 0;
00886 p[3] = 0;
00887 }
00888 return stamp;
00889 }
00890
00891 void MrcIO::swap_header(MrcHeader& mrch)
00892 {
00893 ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
00894 ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
00895 }