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 isFEI(false), is_ri(0), is_new_file(false), initialized(false),
00054 is_transpose(false)
00055 {
00056 memset(&mrch, 0, sizeof(MrcHeader));
00057 is_big_endian = ByteOrder::is_host_big_endian();
00058 }
00059
00060 MrcIO::~MrcIO()
00061 {
00062 if (mrcfile) {
00063 fclose(mrcfile);
00064 mrcfile = 0;
00065 }
00066 }
00067
00068 void MrcIO::init()
00069 {
00070 ENTERFUNC;
00071
00072 if (initialized) {
00073 return;
00074 }
00075
00076 initialized = true;
00077 mrcfile = sfopen(filename, rw_mode, &is_new_file);
00078
00079 if (!is_new_file) {
00080 if (fread(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00081 throw ImageReadException(filename, "MRC header");
00082 }
00083
00084 if (!is_valid(&mrch)) {
00085 throw ImageReadException(filename, "invalid MRC");
00086 }
00087
00088 is_big_endian = ByteOrder::is_data_big_endian(&mrch.nz);
00089 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00090 swap_header(mrch);
00091 }
00092
00093
00094 mode_size = get_mode_size(mrch.mode);
00095 if(is_complex_mode()) {
00096 is_ri = 1;
00097 }
00098
00099 if (mrch.nxstart != 0 || mrch.nystart != 0 || mrch.nzstart != 0) {
00100 LOGWARN("nx/ny/nz start not zero");
00101 }
00102
00103 if (is_complex_mode()) {
00104 mrch.nx *= 2;
00105 }
00106
00107 if (mrch.xlen == 0) {
00108 mrch.xlen = 1.0;
00109 }
00110
00111 if (mrch.ylen == 0) {
00112 mrch.ylen = 1.0;
00113 }
00114
00115 if (mrch.zlen == 0) {
00116 mrch.zlen = 1.0;
00117 }
00118
00119 if(mrch.nlabels>0) {
00120 if( string(mrch.labels[0],3) == "Fei") {
00121 isFEI = true;
00122 }
00123 }
00124
00125 if(mrch.mapc==2 && mrch.mapr==1) {
00126 is_transpose = true;
00127 }
00128 }
00129 EXITFUNC;
00130 }
00131
00132
00133 bool MrcIO::is_image_big_endian()
00134 {
00135 init();
00136 return is_big_endian;
00137 }
00138
00139 bool MrcIO::is_valid(const void *first_block, off_t file_size)
00140 {
00141 ENTERFUNC;
00142
00143 if (!first_block) {
00144 return false;
00145 }
00146
00147 const int *data = static_cast < const int *>(first_block);
00148 int nx = data[0];
00149 int ny = data[1];
00150 int nz = data[2];
00151 int mrcmode = data[3];
00152 int nsymbt = data[23];
00153
00154 bool data_big_endian = ByteOrder::is_data_big_endian(&nz);
00155
00156 if (data_big_endian != ByteOrder::is_host_big_endian()) {
00157 ByteOrder::swap_bytes(&nx);
00158 ByteOrder::swap_bytes(&ny);
00159 ByteOrder::swap_bytes(&nz);
00160 ByteOrder::swap_bytes(&mrcmode);
00161 ByteOrder::swap_bytes(&nsymbt);
00162 }
00163
00164 if (mrcmode == MRC_SHORT_COMPLEX || mrcmode == MRC_FLOAT_COMPLEX) {
00165 nx *= 2;
00166 }
00167
00168 const int max_dim = 1 << 20;
00169
00170 if ((mrcmode >= MRC_UCHAR && mrcmode < MRC_UNKNOWN) &&
00171 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim) && (nz > 0 && nz < max_dim)) {
00172
00173 if (file_size > 0) {
00174 off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz * (off_t)get_mode_size(mrcmode) + (off_t)sizeof(MrcHeader) + nsymbt;
00175 if (file_size == file_size1) {
00176 return true;
00177 }
00178
00179 LOGWARN("image size check fails, still try to read it...");
00180 }
00181 else {
00182 return true;
00183 }
00184
00185 return true;
00186 }
00187 EXITFUNC;
00188 return false;
00189 }
00190
00191 int MrcIO::read_header(Dict & dict, int image_index, const Region * area, bool is_3d)
00192 {
00193 init();
00194
00195 if(isFEI) {
00196 return read_fei_header(dict, image_index, area, is_3d);
00197 }
00198 else {
00199 return read_mrc_header(dict, image_index, area, is_3d);
00200 }
00201 }
00202
00203 int MrcIO::read_mrc_header(Dict & dict, int image_index, const Region * area, bool)
00204 {
00205 ENTERFUNC;
00206
00207
00208 if(image_index < 0) {
00209 image_index = 0;
00210 }
00211 if(image_index != 0) {
00212 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().");
00213 }
00214
00215 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00216
00217 int xlen = 0, ylen = 0, zlen = 0;
00218 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00219
00220 dict["nx"] = xlen;
00221 dict["ny"] = ylen;
00222 dict["nz"] = zlen;
00223 dict["MRC.nx"] = mrch.nx;
00224 dict["MRC.ny"] = mrch.ny;
00225 dict["MRC.nz"] = mrch.nz;
00226
00227 dict["datatype"] = to_em_datatype(mrch.mode);
00228
00229 dict["MRC.nxstart"] = mrch.nxstart;
00230 dict["MRC.nystart"] = mrch.nystart;
00231 dict["MRC.nzstart"] = mrch.nzstart;
00232
00233 dict["MRC.mx"] = mrch.mx;
00234 dict["MRC.my"] = mrch.my;
00235 dict["MRC.mz"] = mrch.mz;
00236
00237 dict["MRC.xlen"] = mrch.xlen;
00238 dict["MRC.ylen"] = mrch.ylen;
00239 dict["MRC.zlen"] = mrch.zlen;
00240
00241 dict["MRC.alpha"] = mrch.alpha;
00242 dict["MRC.beta"] = mrch.beta;
00243 dict["MRC.gamma"] = mrch.gamma;
00244
00245 dict["MRC.mapc"] = mrch.mapc;
00246 dict["MRC.mapr"] = mrch.mapr;
00247 dict["MRC.maps"] = mrch.maps;
00248
00249 dict["MRC.minimum"] = mrch.amin;
00250 dict["MRC.maximum"] = mrch.amax;
00251 dict["MRC.mean"] = mrch.amean;
00252 dict["mean"] = mrch.amean;
00253
00254 dict["MRC.ispg"] = mrch.ispg;
00255 dict["MRC.nsymbt"] = mrch.nsymbt;
00256
00257 float apx = mrch.xlen / mrch.mx;
00258 float apy = mrch.ylen / mrch.my;
00259 float apz = mrch.zlen / mrch.mz;
00260 if(apx>1000 || apx<0.01) {
00261 dict["apix_x"] = 1.0f;
00262 }
00263 else {
00264 dict["apix_x"] = apx;
00265 }
00266 if(apy>1000 || apy<0.01) {
00267 dict["apix_y"] = 1.0f;
00268 }
00269 else {
00270 dict["apix_y"] = apy;
00271 }
00272 if(apz>1000 || apz<0.01) {
00273 dict["apix_z"] = 1.0f;
00274 }
00275 else {
00276 dict["apix_z"] = apz;
00277 }
00278
00279 if (area) {
00280 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
00281 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
00282
00283 if (area->get_ndim() == 3 && mrch.nz > 1) {
00284 dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
00285 }
00286 else {
00287 dict["origin_z"] = mrch.zorigin;
00288 }
00289 }
00290 else {
00291 dict["origin_x"] = mrch.xorigin;
00292 dict["origin_y"] = mrch.yorigin;
00293 dict["origin_z"] = mrch.zorigin;
00294 }
00295
00296 if (is_complex_mode()) {
00297 dict["is_complex"] = 1;
00298 dict["is_complex_ri"] = 1;
00299 }
00300
00301 dict["MRC.machinestamp"] = mrch.machinestamp;
00302
00303 dict["MRC.rms"] = mrch.rms;
00304 dict["sigma"] = mrch.rms;
00305 dict["MRC.nlabels"] = mrch.nlabels;
00306 for (int i = 0; i < mrch.nlabels; i++) {
00307 char label[32];
00308 sprintf(label, "MRC.label%d", i);
00309 dict[string(label)] = string(mrch.labels[i],80);
00310 }
00311
00312 EMAN1Ctf ctf_;
00313 if(read_ctf(ctf_) == 0) {
00314 vector<float> vctf = ctf_.to_vector();
00315 dict["ctf"] = vctf;
00316 }
00317
00318 if(is_transpose) {
00319 dict["nx"] = ylen;
00320 dict["ny"] = xlen;
00321 dict["MRC.nx"] = mrch.ny;
00322 dict["MRC.ny"] = mrch.nx;
00323 dict["MRC.mx"] = mrch.my;
00324 dict["MRC.my"] = mrch.mx;
00325 dict["apix_x"] = mrch.ylen / mrch.my;
00326 dict["apix_y"] = mrch.xlen / mrch.mx;
00327 dict["origin_x"] = mrch.yorigin;
00328 dict["origin_y"] = mrch.xorigin;
00329 dict["MRC.nxstart"] = mrch.nystart;
00330 dict["MRC.nystart"] = mrch.nxstart;
00331 }
00332
00333 Transform * trans = new Transform();
00334 if(is_transpose) {
00335 trans->set_trans(mrch.nystart, mrch.nxstart, mrch.nzstart);
00336 }
00337 else {
00338 trans->set_trans(mrch.nxstart, mrch.nystart, mrch.nzstart);
00339 }
00340
00341 if(zlen<=1) {
00342 dict["xform.projection"] = trans;
00343 }
00344 else {
00345 dict["xform.align3d"] = trans;
00346 }
00347
00348 if(trans) {delete trans; trans=0;}
00349
00350 EXITFUNC;
00351 return 0;
00352 }
00353
00354 int MrcIO::read_fei_header(Dict & dict, int image_index, const Region * area, bool)
00355 {
00356 ENTERFUNC;
00357
00358 if(image_index < 0) {
00359 image_index = 0;
00360 }
00361
00362 if(area && area->get_depth() > 1) {
00363 throw ImageDimensionException("FEI MRC image is 2D tile series, only 2D regional reading accepted.");
00364 }
00365
00366 init();
00367
00368 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file,false);
00369
00370 int xlen = 0, ylen = 0, zlen = 0;
00371 EMUtil::get_region_dims(area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
00372
00373 dict["nx"] = xlen;
00374 dict["ny"] = ylen;
00375 dict["nz"] = 1;
00376 dict["FEIMRC.nx"] = feimrch.nx;
00377 dict["FEIMRC.ny"] = feimrch.ny;
00378 dict["FEIMRC.nz"] = feimrch.nz;
00379
00380 dict["datatype"] = to_em_datatype(feimrch.mode);
00381
00382 dict["FEIMRC.nxstart"] = feimrch.nxstart;
00383 dict["FEIMRC.nystart"] = feimrch.nystart;
00384 dict["FEIMRC.nzstart"] = feimrch.nzstart;
00385
00386 dict["FEIMRC.mx"] = feimrch.mx;
00387 dict["FEIMRC.my"] = feimrch.my;
00388 dict["FEIMRC.mz"] = feimrch.mz;
00389
00390 dict["FEIMRC.xlen"] = feimrch.xlen;
00391 dict["FEIMRC.ylen"] = feimrch.ylen;
00392 dict["FEIMRC.zlen"] = feimrch.zlen;
00393
00394 dict["FEIMRC.alpha"] = feimrch.alpha;
00395 dict["FEIMRC.beta"] = feimrch.beta;
00396 dict["FEIMRC.gamma"] = feimrch.gamma;
00397
00398 dict["FEIMRC.mapc"] = feimrch.mapc;
00399 dict["FEIMRC.mapr"] = feimrch.mapr;
00400 dict["FEIMRC.maps"] = feimrch.maps;
00401
00402 dict["FEIMRC.minimum"] = feimrch.amin;
00403 dict["FEIMRC.maximum"] = feimrch.amax;
00404 dict["FEIMRC.mean"] = feimrch.amean;
00405 dict["mean"] = feimrch.amean;
00406
00407 dict["FEIMRC.ispg"] = feimrch.ispg;
00408 dict["FEIMRC.nsymbt"] = feimrch.nsymbt;
00409
00410 dict["apix_x"] = feimrch.xlen / feimrch.mx;
00411 dict["apix_y"] = feimrch.ylen / feimrch.my;
00412 dict["apix_z"] = feimrch.zlen / feimrch.mz;
00413
00414 dict["FEIMRC.next"] = feimrch.next;
00415 dict["FEIMRC.dvid"] = feimrch.dvid;
00416 dict["FEIMRC.numintegers"] = feimrch.numintegers;
00417 dict["FEIMRC.numfloats"] = feimrch.numfloats;
00418 dict["FEIMRC.sub"] = feimrch.sub;
00419 dict["FEIMRC.zfac"] = feimrch.zfac;
00420
00421 dict["FEIMRC.min2"] = feimrch.min2;
00422 dict["FEIMRC.max2"] = feimrch.max2;
00423 dict["FEIMRC.min3"] = feimrch.min3;
00424 dict["FEIMRC.max3"] = feimrch.max3;
00425 dict["FEIMRC.min4"] = feimrch.min4;
00426 dict["FEIMRC.max4"] = feimrch.max4;
00427
00428 dict["FEIMRC.idtype"] = feimrch.idtype;
00429 dict["FEIMRC.lens"] = feimrch.lens;
00430 dict["FEIMRC.nd1"] = feimrch.nd1;
00431 dict["FEIMRC.nd2"] = feimrch.nd2;
00432 dict["FEIMRC.vd1"] = feimrch.vd1;
00433 dict["FEIMRC.vd2"] = feimrch.vd2;
00434
00435 for(int i=0; i<9; i++) {
00436 char label[32];
00437 sprintf(label, "MRC.tiltangles%d", i);
00438 dict[string(label)] = feimrch.tiltangles[i];
00439 }
00440
00441 dict["FEIMRC.zorg"] = feimrch.zorg;
00442 dict["FEIMRC.xorg"] = feimrch.xorg;
00443 dict["FEIMRC.yorg"] = feimrch.yorg;
00444
00445 dict["FEIMRC.nlabl"] = feimrch.nlabl;
00446 for (int i = 0; i < feimrch.nlabl; i++) {
00447 char label[32];
00448 sprintf(label, "MRC.label%d", i);
00449 dict[string(label)] = string(feimrch.labl[i], 80);
00450 }
00451
00452
00453 FeiMrcExtHeader feiexth;
00454 portable_fseek(mrcfile, sizeof(FeiMrcHeader)+sizeof(FeiMrcExtHeader)*image_index, SEEK_SET);
00455 if (fread(&feiexth, sizeof(FeiMrcExtHeader), 1, mrcfile) != 1) {
00456 throw ImageReadException(filename, "FEI MRC extended header");
00457 }
00458
00459 dict["FEIMRC.a_tilt"] = feiexth.a_tilt;
00460 dict["FEIMRC.b_tilt"] = feiexth.b_tilt;
00461
00462 dict["FEIMRC.x_stage"] = feiexth.x_stage;
00463 dict["FEIMRC.y_stage"] = feiexth.y_stage;
00464 dict["FEIMRC.z_stage"] = feiexth.z_stage;
00465
00466 dict["FEIMRC.x_shift"] = feiexth.x_shift;
00467 dict["FEIMRC.y_shift"] = feiexth.y_shift;
00468
00469 dict["FEIMRC.defocus"] = feiexth.defocus;
00470 dict["FEIMRC.exp_time"] = feiexth.exp_time;
00471 dict["FEIMRC.mean_int"] = feiexth.mean_int;
00472 dict["FEIMRC.tilt_axis"] = feiexth.tilt_axis;
00473
00474 dict["FEIMRC.pixel_size"] = feiexth.pixel_size;
00475 dict["FEIMRC.magnification"] = feiexth.magnification;
00476 dict["FEIMRC.ht"] = feiexth.ht;
00477 dict["FEIMRC.binning"] = feiexth.binning;
00478 dict["FEIMRC.appliedDefocus"] = feiexth.appliedDefocus;
00479
00480
00481
00482 EXITFUNC;
00483 return 0;
00484 }
00485
00486 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00487 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00488 {
00489 ENTERFUNC;
00490
00491
00492 if(image_index == -1) {
00493 image_index = 0;
00494 }
00495 if(image_index != 0) {
00496 throw ImageWriteException(filename, "MRC file does not support stack.");
00497 }
00498 check_write_access(rw_mode, image_index, 1);
00499 if (area) {
00500 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00501 EXITFUNC;
00502 return 0;
00503 }
00504
00505 int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00506 int nx = dict["nx"];
00507 int ny = dict["ny"];
00508 int nz = dict["nz"];
00509 is_ri = dict["is_complex_ri"];
00510
00511 bool opposite_endian = false;
00512
00513 if (!is_new_file) {
00514 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00515 opposite_endian = true;
00516 }
00517 #if 0
00518 if (new_mode != mrch.mode) {
00519 LOGERR("cannot write to different mode file %s", filename.c_str());
00520 return 1;
00521 }
00522 #endif
00523 portable_fseek(mrcfile, 0, SEEK_SET);
00524 }
00525 else {
00526 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00527 mrch.mapc = 1;
00528 mrch.mapr = 2;
00529 mrch.maps = 3;
00530 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00531 }
00532
00533 if(nz<=1 && dict.has_key("xform.projection") && !dict.has_key("UCSF.chimera")) {
00534 Transform * t = dict["xform.projection"];
00535 Dict d = t->get_params("imagic");
00536 mrch.alpha = d["alpha"];
00537 mrch.beta = d["beta"];
00538 mrch.gamma = d["gamma"];
00539 mrch.xorigin = d["tx"];
00540 mrch.yorigin = d["ty"];
00541 mrch.zorigin = d["tz"];
00542 if(t) {delete t; t=0;}
00543 }
00544 else if(nz>1 && dict.has_key("xform.align3d") && !dict.has_key("UCSF.chimera")) {
00545 Transform * t = dict["xform.align3d"];
00546 Dict d = t->get_params("imagic");
00547 mrch.alpha = d["alpha"];
00548 mrch.beta = d["beta"];
00549 mrch.gamma = d["gamma"];
00550 mrch.xorigin = d["tx"];
00551 mrch.yorigin = d["ty"];
00552 mrch.zorigin = d["tz"];
00553 if(t) {delete t; t=0;}
00554 }
00555
00556 if(dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00557 mrch.xorigin = (float)dict["origin_x"];
00558 mrch.yorigin = (float)dict["origin_y"];
00559
00560 if (is_new_file) {
00561 mrch.zorigin = (float)dict["origin_z"];
00562 }
00563 else {
00564 mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00565 }
00566 }
00567
00568 if (dict.has_key("MRC.nlabels")) {
00569 mrch.nlabels = dict["MRC.nlabels"];
00570 }
00571
00572 for (int i = 0; i < MRC_NUM_LABELS; i++) {
00573 char label[32];
00574 #ifdef _WIN32
00575 _snprintf(label,31, "MRC.label%d", i);
00576 #else
00577 snprintf(label,31, "MRC.label%d", i);
00578 #endif //_WIN32
00579 if (dict.has_key(label)) {
00580 #ifdef _WIN32
00581 _snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00582 #else
00583 snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00584 #endif //_WIN32
00585 mrch.nlabels = i + 1;
00586 }
00587 }
00588
00589 if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00590 #ifdef _WIN32
00591 _snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00592 #else
00593 snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00594 #endif //_WIN32
00595 mrch.nlabels++;
00596 }
00597
00598 mrch.labels[mrch.nlabels][0] = '\0';
00599 mrch.mode = new_mode;
00600
00601 if (is_complex_mode()) {
00602 mrch.nx = nx / 2;
00603 }
00604 else {
00605 mrch.nx = nx;
00606 }
00607 mrch.ny = ny;
00608
00609 if (is_new_file) {
00610 mrch.nz = nz;
00611 }
00612 else if (image_index >= mrch.nz) {
00613 mrch.nz = image_index + 1;
00614 }
00615
00616 mrch.ispg = dict.has_key("MRC.ispg") ? (int)dict["MRC.ispg"] : 0;
00617 mrch.nsymbt = 0;
00618 mrch.amin = dict["minimum"];
00619 mrch.amax = dict["maximum"];
00620 mrch.amean = dict["mean"];
00621 mrch.rms = dict["sigma"];
00622
00625
00626
00627
00628
00629 mrch.mx = nx;
00630
00631
00632
00633
00634
00635 mrch.my = ny;
00636
00637
00638
00639
00640
00641 mrch.mz = nz;
00642
00643
00644 mrch.xlen = mrch.mx * (float) dict["apix_x"];
00645 mrch.ylen = mrch.my * (float) dict["apix_y"];
00646 mrch.zlen = mrch.mz * (float) dict["apix_z"];
00647
00648 if(dict.has_key("MRC.nxstart")) {
00649 mrch.nxstart = dict["MRC.nxstart"];
00650 }
00651 else {
00652 mrch.nxstart = -nx / 2;
00653 }
00654 if(dict.has_key("MRC.nystart")) {
00655 mrch.nystart = dict["MRC.nystart"];
00656 }
00657 else {
00658 mrch.nystart = -ny / 2;
00659 }
00660 if(dict.has_key("MRC.nzstart")) {
00661 mrch.nzstart = dict["MRC.nzstart"];
00662 }
00663 else {
00664 mrch.nzstart = -nz / 2;
00665 }
00666
00667 strncpy(mrch.map,"MAP ",4);
00668 mrch.machinestamp = generate_machine_stamp();
00669
00670 MrcHeader mrch2 = mrch;
00671
00672 if (opposite_endian || !use_host_endian) {
00673 swap_header(mrch2);
00674 }
00675
00676 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00677 throw ImageWriteException(filename, "MRC header");
00678 }
00679
00680 mode_size = get_mode_size(mrch.mode);
00681 is_new_file = false;
00682
00683 if( dict.has_key("ctf") ) {
00684 vector<float> vctf = dict["ctf"];
00685 EMAN1Ctf ctf_;
00686 ctf_.from_vector(vctf);
00687 write_ctf(ctf_);
00688 }
00689
00690 EXITFUNC;
00691 return 0;
00692 }
00693
00694 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool )
00695 {
00696 ENTERFUNC;
00697
00698 if(!isFEI) {
00699
00700 image_index = 0;
00701 }
00702
00703 if(is_transpose && area!=0) {
00704 printf("Warning: This image dimension is in (y,x,z), region I/O not supported, return the whole image instead.");
00705 }
00706
00707 check_read_access(image_index, rdata);
00708
00709 if (area && is_complex_mode()) {
00710 LOGERR("Error: cannot read a region of a complex image.");
00711 return 1;
00712 }
00713
00714 unsigned char *cdata = (unsigned char *) rdata;
00715 short *sdata = (short *) rdata;
00716 unsigned short *usdata = (unsigned short *) rdata;
00717
00718 size_t size = 0;
00719 int xlen = 0, ylen = 0, zlen = 0;
00720 if(isFEI) {
00721 check_region(area, FloatSize(mrch.nx, mrch.ny, 1), is_new_file, false);
00722 portable_fseek(mrcfile, sizeof(MrcHeader)+feimrch.next, SEEK_SET);
00723
00724 Region * new_area;
00725 if(area) {
00726 new_area = new Region(area->x_origin(), area->y_origin(), (float)image_index, area->get_width(), area->get_height(), 1.0f);
00727 }
00728 else {
00729 new_area = new Region(0, 0, image_index, feimrch.nx, feimrch.ny, 1);
00730 }
00731 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00732 image_index, mode_size,
00733 feimrch.nx, feimrch.ny, feimrch.nz, new_area);
00734
00735 EMUtil::get_region_dims(new_area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
00736
00737 size = (size_t)xlen * ylen * zlen;
00738
00739 delete new_area;
00740 }
00741 else {
00742 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00743 portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00744
00745 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00746 image_index, mode_size,
00747 mrch.nx, mrch.ny, mrch.nz, area);
00748
00749 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00750
00751 size = xlen * ylen * zlen;
00752 }
00753
00754 if (mrch.mode != MRC_UCHAR) {
00755 if (mode_size == sizeof(short)) {
00756 become_host_endian < short >(sdata, size);
00757 }
00758 else if (mode_size == sizeof(float)) {
00759 become_host_endian < float >(rdata, size);
00760 }
00761 }
00762
00763 if (mrch.mode == MRC_UCHAR) {
00764 for (size_t i = 0; i < size; ++i) {
00765 size_t j = size - 1 - i;
00766
00767 rdata[j] = static_cast < float >(cdata[j]);
00768 }
00769 }
00770 else if (mrch.mode == MRC_SHORT ) {
00771 for (size_t i = 0; i < size; ++i) {
00772 size_t j = size - 1 - i;
00773 rdata[j] = static_cast < float >(sdata[j]);
00774 }
00775 }
00776 else if (mrch.mode == MRC_USHORT) {
00777 for (size_t i = 0; i < size; ++i) {
00778 size_t j = size - 1 - i;
00779 rdata[j] = static_cast < float >(usdata[j]);
00780 }
00781 }
00782
00783 if(is_transpose) {
00784 transpose(rdata, xlen, ylen, zlen);
00785 }
00786
00787 if (is_complex_mode()) {
00788 if(!is_ri) Util::ap2ri(rdata, size);
00789 Util::flip_complex_phase(rdata, size);
00790 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00791 }
00792 EXITFUNC;
00793 return 0;
00794 }
00795
00796 int MrcIO::write_data(float *data, int image_index, const Region* area,
00797 EMUtil::EMDataType, bool use_host_endian)
00798 {
00799 ENTERFUNC;
00800
00801 image_index = 0;
00802 check_write_access(rw_mode, image_index, 1, data);
00803 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00804
00805 int nx = mrch.nx;
00806 int ny = mrch.ny;
00807 int nz = mrch.nz;
00808 size_t size = (size_t)nx * ny * nz;
00809
00810 if (is_complex_mode()) {
00811 nx *= 2;
00812 if (!is_ri) {
00813 Util::ap2ri(data, size);
00814 is_ri = 1;
00815 }
00816 Util::flip_complex_phase(data, size);
00817 Util::rotate_phase_origin(data, nx, ny, nz);
00818 }
00819
00820 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00821
00822 if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian) {
00823 if (mrch.mode != MRC_UCHAR) {
00824 if (mode_size == sizeof(short)) {
00825 ByteOrder::swap_bytes((short*) data, size);
00826 }
00827 else if (mode_size == sizeof(float)) {
00828 ByteOrder::swap_bytes((float*) data, size);
00829 }
00830 }
00831 }
00832 mode_size = get_mode_size(mrch.mode);
00833
00834
00835
00836
00837 void * ptr_data = data;
00838
00839 float rendermin = 0.0f;
00840 float rendermax = 0.0f;
00841 EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, nz);
00842
00843 unsigned char *cdata = 0;
00844 short *sdata = 0;
00845 unsigned short *usdata = 0;
00846 if (mrch.mode == MRC_UCHAR) {
00847 cdata = new unsigned char[size];
00848 for (size_t i = 0; i < size; ++i) {
00849 if(data[i] <= rendermin) {
00850 cdata[i] = 0;
00851 }
00852 else if(data[i] >= rendermax){
00853 cdata[i] = UCHAR_MAX;
00854 }
00855 else {
00856 cdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
00857 }
00858 }
00859 ptr_data = cdata;
00860 update_stat((void *)cdata);
00861 }
00862 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00863 sdata = new short[size];
00864 for (size_t i = 0; i < size; ++i) {
00865 if(data[i] <= rendermin) {
00866 sdata[i] = SHRT_MIN;
00867 }
00868 else if(data[i] >= rendermax) {
00869 sdata[i] = SHRT_MAX;
00870 }
00871 else {
00872 sdata[i]=(short)(((data[i]-rendermin)/(rendermax-rendermin))*(SHRT_MAX-SHRT_MIN) - SHRT_MAX);
00873 }
00874 }
00875 ptr_data = sdata;
00876 update_stat((void *)sdata);
00877 }
00878 else if (mrch.mode == MRC_USHORT) {
00879 usdata = new unsigned short[size];
00880 for (size_t i = 0; i < size; ++i) {
00881 if(data[i] <= rendermin) {
00882 usdata[i] = 0;
00883 }
00884 else if(data[i] >= rendermax) {
00885 usdata[i] = USHRT_MAX;
00886 }
00887 else {
00888 usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
00889 }
00890 }
00891 ptr_data = usdata;
00892 update_stat((void *)usdata);
00893 }
00894
00895
00896
00897
00898 EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
00899 mode_size, nx, mrch.ny, mrch.nz, area);
00900
00901 if(cdata) {delete [] cdata; cdata=0;}
00902 if(sdata) {delete [] sdata; sdata=0;}
00903 if(usdata) {delete [] usdata; usdata=0;}
00904
00905 #if 0
00906 int row_size = nx * get_mode_size(mrch.mode);
00907 int sec_size = nx * ny;
00908
00909 unsigned char *cbuf = new unsigned char[row_size];
00910 unsigned short *sbuf = (unsigned short *) cbuf;
00911
00912 for (int i = 0; i < nz; i++) {
00913 int i2 = i * sec_size;
00914 for (int j = 0; j < ny; j++) {
00915 int k = i2 + j * nx;
00916 void *pbuf = 0;
00917
00918 switch (mrch.mode) {
00919 case MRC_UCHAR:
00920 for (int l = 0; l < nx; l++) {
00921 cbuf[l] = static_cast < unsigned char >(data[k + l]);
00922 }
00923 pbuf = cbuf;
00924 fwrite(cbuf, row_size, 1, mrcfile);
00925 break;
00926
00927 case MRC_SHORT:
00928 case MRC_SHORT_COMPLEX:
00929 for (int l = 0; l < nx; l++) {
00930 sbuf[l] = static_cast < short >(data[k + l]);
00931 }
00932 pbuf = sbuf;
00933 fwrite(sbuf, row_size, 1, mrcfile);
00934 break;
00935
00936 case MRC_USHORT:
00937 for (int l = 0; l < nx; l++) {
00938 sbuf[l] = static_cast < unsigned short >(data[k + l]);
00939 }
00940 pbuf = sbuf;
00941 fwrite(sbuf, row_size, 1, mrcfile);
00942 break;
00943
00944 case MRC_FLOAT:
00945 case MRC_FLOAT_COMPLEX:
00946 pbuf = &data[k];
00947 break;
00948 }
00949 if (pbuf) {
00950 fwrite(pbuf, row_size, 1, mrcfile);
00951 }
00952 }
00953 }
00954
00955 if(cbuf)
00956 {
00957 delete[]cbuf;
00958 cbuf = 0;
00959 }
00960 #endif
00961
00962 EXITFUNC;
00963 return 0;
00964 }
00965
00966 void MrcIO::update_stat(void* data)
00967 {
00968 size_t size = mrch.nx * mrch.ny * mrch.nz;
00969 float v = 0.0f;
00970 double sum = 0.0;
00971 double square_sum = 0.0;
00972 double mean = 0.0;
00973 float min, max;
00974
00975 unsigned char * cdata = 0;
00976 short * sdata = 0;
00977 unsigned short * usdata = 0;
00978
00979 if (mrch.mode == MRC_UCHAR) {
00980 max = 0.0f;
00981 min = UCHAR_MAX;
00982 cdata = (unsigned char *)data;
00983
00984 for (size_t i = 0; i < size; ++i) {
00985 v = (float)(cdata[i]);
00986 #ifdef _WIN32
00987 max = _cpp_max(max,v);
00988 min = _cpp_min(min,v);
00989 #else
00990 max=std::max<float>(max,v);
00991 min=std::min<float>(min,v);
00992 #endif //_WIN32
00993
00994 sum += v;
00995 square_sum += v * v;
00996 }
00997 }
00998 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00999 max = (float)SHRT_MIN;
01000 min = (float)SHRT_MAX;
01001 sdata = (short *)data;
01002
01003 for (size_t i = 0; i < size; ++i) {
01004 v = (float)(sdata[i]);
01005 #ifdef _WIN32
01006 max = _cpp_max(max,v);
01007 min = _cpp_min(min,v);
01008 #else
01009 max=std::max<float>(max,v);
01010 min=std::min<float>(min,v);
01011 #endif //_WIN32
01012
01013 sum += v;
01014 square_sum += v * v;
01015 }
01016 }
01017 else if (mrch.mode == MRC_USHORT) {
01018 max = 0.0f;
01019 min = (float)USHRT_MAX;
01020 usdata = (unsigned short*)data;
01021
01022 for (size_t i = 0; i < size; ++i) {
01023 v = (float)(usdata[i]);
01024 #ifdef _WIN32
01025 max = _cpp_max(max,v);
01026 min = _cpp_min(min,v);
01027 #else
01028 max=std::max<float>(max,v);
01029 min=std::min<float>(min,v);
01030 #endif //_WIN32
01031
01032 sum += v;
01033 square_sum += v * v;
01034 }
01035 }
01036 else {
01037 throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
01038 }
01039
01040 mean = sum/size;
01041 #ifdef _WIN32
01042 float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / size)/(size-1)));
01043 #else
01044 float sigma = (float)std::sqrt(std::max<float>(0.0,(square_sum - sum*sum / size)/(size-1)));
01045 #endif //_WIN32
01046
01047
01048 mrch.amin = min;
01049 mrch.amax = max;
01050 mrch.amean = (float)mean;
01051 mrch.rms = sigma;
01052
01053 MrcHeader mrch2 = mrch;
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068
01069
01070 portable_fseek(mrcfile, 0, SEEK_SET);
01071
01072 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
01073 throw ImageWriteException(filename, "MRC header");
01074 }
01075
01076 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
01077 }
01078
01079 bool MrcIO::is_complex_mode()
01080 {
01081 init();
01082 if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
01083 return true;
01084 }
01085 return false;
01086 }
01087
01088
01089 int MrcIO::read_ctf(Ctf & ctf, int)
01090 {
01091 ENTERFUNC;
01092 init();
01093 size_t n = strlen(CTF_MAGIC);
01094
01095 int err = 1;
01096 if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
01097 err = ctf.from_string(string(&mrch.labels[0][n]));
01098 }
01099 EXITFUNC;
01100 return err;
01101 }
01102
01103 void MrcIO::write_ctf(const Ctf & ctf, int)
01104 {
01105 ENTERFUNC;
01106 init();
01107
01108 string ctf_str = ctf.to_string();
01109 #ifdef _WIN32
01110 _snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01111 #else
01112 snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01113 #endif //_WIN32
01114 rewind(mrcfile);
01115
01116 if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
01117 throw ImageWriteException(filename, "write CTF info to header failed");
01118 }
01119 EXITFUNC;
01120 }
01121
01122 void MrcIO::flush()
01123 {
01124 fflush(mrcfile);
01125 }
01126
01127
01128 int MrcIO::get_mode_size(int mm)
01129 {
01130 MrcIO::MrcMode m = static_cast < MrcMode > (mm);
01131
01132 int msize = 0;
01133 switch (m) {
01134 case MRC_UCHAR:
01135 msize = sizeof(char);
01136 break;
01137 case MRC_SHORT:
01138 case MRC_USHORT:
01139 case MRC_SHORT_COMPLEX:
01140 msize = sizeof(short);
01141 break;
01142 case MRC_FLOAT:
01143 case MRC_FLOAT_COMPLEX:
01144 msize = sizeof(float);
01145 break;
01146 default:
01147 msize = 0;
01148 }
01149
01150 return msize;
01151 }
01152
01153 int MrcIO::to_em_datatype(int m)
01154 {
01155 EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
01156
01157 switch (m) {
01158 case MRC_UCHAR:
01159 e = EMUtil::EM_UCHAR;
01160 break;
01161 case MRC_SHORT:
01162 e = EMUtil::EM_SHORT;
01163 break;
01164 case MRC_USHORT:
01165 e = EMUtil::EM_USHORT;
01166 break;
01167 case MRC_SHORT_COMPLEX:
01168 e = EMUtil::EM_SHORT_COMPLEX;
01169 break;
01170 case MRC_FLOAT:
01171 e = EMUtil::EM_FLOAT;
01172 break;
01173 case MRC_FLOAT_COMPLEX:
01174 e = EMUtil::EM_FLOAT_COMPLEX;
01175 break;
01176 default:
01177 e = EMUtil::EM_UNKNOWN;
01178 }
01179 return e;
01180 }
01181
01182
01183 int MrcIO::to_mrcmode(int e, int is_complex)
01184 {
01185 MrcMode m = MRC_UNKNOWN;
01186 EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
01187
01188 switch (em_type) {
01189 case EMUtil::EM_UCHAR:
01190 m = MRC_UCHAR;
01191 break;
01192 case EMUtil::EM_USHORT:
01193 if (is_complex) {
01194 m = MRC_SHORT_COMPLEX;
01195 }
01196 else {
01197 m = MRC_USHORT;
01198 }
01199 break;
01200 case EMUtil::EM_SHORT:
01201 if (is_complex) {
01202 m = MRC_SHORT_COMPLEX;
01203 }
01204 else {
01205 m = MRC_SHORT;
01206 }
01207 break;
01208 case EMUtil::EM_SHORT_COMPLEX:
01209 case EMUtil::EM_USHORT_COMPLEX:
01210 m = MRC_SHORT_COMPLEX;
01211 break;
01212 case EMUtil::EM_CHAR:
01213 case EMUtil::EM_INT:
01214 case EMUtil::EM_UINT:
01215 case EMUtil::EM_FLOAT:
01216 if (is_complex) {
01217 m = MRC_FLOAT_COMPLEX;
01218 }
01219 else {
01220 m = MRC_FLOAT;
01221 }
01222 break;
01223 case EMUtil::EM_FLOAT_COMPLEX:
01224 m = MRC_FLOAT_COMPLEX;
01225 break;
01226 default:
01227 m = MRC_FLOAT;
01228 }
01229
01230 return m;
01231 }
01232
01233
01234
01235 int MrcIO::generate_machine_stamp()
01236 {
01237 int stamp = 0;
01238 char *p = (char *) (&stamp);
01239
01240 if (ByteOrder::is_host_big_endian()) {
01241 p[0] = 0x11;
01242 p[1] = 0x11;
01243 p[2] = 0;
01244 p[3] = 0;
01245 }
01246 else {
01247 p[0] = 0x44;
01248 p[1] = 0x41;
01249 p[2] = 0;
01250 p[3] = 0;
01251 }
01252 return stamp;
01253 }
01254
01255 void MrcIO::swap_header(MrcHeader& mrch)
01256 {
01257 ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
01258 ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
01259 }
01260
01261 int MrcIO::get_nimg()
01262 {
01263 init();
01264
01265 if(isFEI) {
01266 return feimrch.nz;
01267 }
01268 else {
01269 return 1;
01270 }
01271 }
01272
01273 int MrcIO::transpose(float *data, int xlen, int ylen, int zlen) const
01274 {
01275 float * tmp = new float[xlen*ylen];
01276
01277 for(size_t z=0; z<(size_t)zlen; ++z) {
01278 for(size_t y=0; y<(size_t)ylen; ++y) {
01279 for(size_t x=0; x<(size_t)xlen; ++x) {
01280 tmp[x*ylen+y] = data[z*xlen*ylen+y*xlen+x];
01281 }
01282 }
01283 std::copy(tmp, tmp+xlen*ylen, data+z*xlen*ylen);
01284 }
01285
01286 delete [] tmp;
01287
01288 return 0;
01289 }