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 trans->set_rotation(Dict("type", "imagic", "alpha", mrch.alpha, "beta", mrch.beta, "gamma", mrch.gamma));
00337 }
00338 else {
00339 trans->set_trans(mrch.nxstart, mrch.nystart, mrch.nzstart);
00340 trans->set_rotation(Dict("type", "imagic", "alpha", mrch.alpha, "beta", mrch.beta, "gamma", mrch.gamma));
00341 }
00342
00343 if(zlen<=1) {
00344 dict["xform.projection"] = trans;
00345 }
00346 else {
00347 dict["xform.align3d"] = trans;
00348 }
00349
00350 if(trans) {delete trans; trans=0;}
00351
00352 EXITFUNC;
00353 return 0;
00354 }
00355
00356 int MrcIO::read_fei_header(Dict & dict, int image_index, const Region * area, bool)
00357 {
00358 ENTERFUNC;
00359
00360 if(image_index < 0) {
00361 image_index = 0;
00362 }
00363
00364 init();
00365
00366 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
00367
00368 int xlen = 0, ylen = 0, zlen = 0;
00369 EMUtil::get_region_dims(area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
00370
00371 dict["nx"] = xlen;
00372 dict["ny"] = ylen;
00373 dict["nz"] = zlen;
00374 dict["FEIMRC.nx"] = feimrch.nx;
00375 dict["FEIMRC.ny"] = feimrch.ny;
00376 dict["FEIMRC.nz"] = feimrch.nz;
00377
00378 dict["datatype"] = to_em_datatype(feimrch.mode);
00379
00380 dict["FEIMRC.nxstart"] = feimrch.nxstart;
00381 dict["FEIMRC.nystart"] = feimrch.nystart;
00382 dict["FEIMRC.nzstart"] = feimrch.nzstart;
00383
00384 dict["FEIMRC.mx"] = feimrch.mx;
00385 dict["FEIMRC.my"] = feimrch.my;
00386 dict["FEIMRC.mz"] = feimrch.mz;
00387
00388 dict["FEIMRC.xlen"] = feimrch.xlen;
00389 dict["FEIMRC.ylen"] = feimrch.ylen;
00390 dict["FEIMRC.zlen"] = feimrch.zlen;
00391
00392 dict["FEIMRC.alpha"] = feimrch.alpha;
00393 dict["FEIMRC.beta"] = feimrch.beta;
00394 dict["FEIMRC.gamma"] = feimrch.gamma;
00395
00396 dict["FEIMRC.mapc"] = feimrch.mapc;
00397 dict["FEIMRC.mapr"] = feimrch.mapr;
00398 dict["FEIMRC.maps"] = feimrch.maps;
00399
00400 dict["FEIMRC.minimum"] = feimrch.amin;
00401 dict["FEIMRC.maximum"] = feimrch.amax;
00402 dict["FEIMRC.mean"] = feimrch.amean;
00403 dict["mean"] = feimrch.amean;
00404
00405 dict["FEIMRC.ispg"] = feimrch.ispg;
00406 dict["FEIMRC.nsymbt"] = feimrch.nsymbt;
00407
00408 dict["apix_x"] = feimrch.xlen / feimrch.mx;
00409 dict["apix_y"] = feimrch.ylen / feimrch.my;
00410 dict["apix_z"] = feimrch.zlen / feimrch.mz;
00411
00412 dict["FEIMRC.next"] = feimrch.next;
00413 dict["FEIMRC.dvid"] = feimrch.dvid;
00414 dict["FEIMRC.numintegers"] = feimrch.numintegers;
00415 dict["FEIMRC.numfloats"] = feimrch.numfloats;
00416 dict["FEIMRC.sub"] = feimrch.sub;
00417 dict["FEIMRC.zfac"] = feimrch.zfac;
00418
00419 dict["FEIMRC.min2"] = feimrch.min2;
00420 dict["FEIMRC.max2"] = feimrch.max2;
00421 dict["FEIMRC.min3"] = feimrch.min3;
00422 dict["FEIMRC.max3"] = feimrch.max3;
00423 dict["FEIMRC.min4"] = feimrch.min4;
00424 dict["FEIMRC.max4"] = feimrch.max4;
00425
00426 dict["FEIMRC.idtype"] = feimrch.idtype;
00427 dict["FEIMRC.lens"] = feimrch.lens;
00428 dict["FEIMRC.nd1"] = feimrch.nd1;
00429 dict["FEIMRC.nd2"] = feimrch.nd2;
00430 dict["FEIMRC.vd1"] = feimrch.vd1;
00431 dict["FEIMRC.vd2"] = feimrch.vd2;
00432
00433 for(int i=0; i<9; i++) {
00434 char label[32];
00435 sprintf(label, "MRC.tiltangles%d", i);
00436 dict[string(label)] = feimrch.tiltangles[i];
00437 }
00438
00439 dict["FEIMRC.zorg"] = feimrch.zorg;
00440 dict["FEIMRC.xorg"] = feimrch.xorg;
00441 dict["FEIMRC.yorg"] = feimrch.yorg;
00442
00443 dict["FEIMRC.nlabl"] = feimrch.nlabl;
00444 for (int i = 0; i < feimrch.nlabl; i++) {
00445 char label[32];
00446 sprintf(label, "MRC.label%d", i);
00447 dict[string(label)] = string(feimrch.labl[i], 80);
00448 }
00449
00450
00451 FeiMrcExtHeader feiexth;
00452 portable_fseek(mrcfile, sizeof(FeiMrcHeader)+sizeof(FeiMrcExtHeader)*image_index, SEEK_SET);
00453 if (fread(&feiexth, sizeof(FeiMrcExtHeader), 1, mrcfile) != 1) {
00454 throw ImageReadException(filename, "FEI MRC extended header");
00455 }
00456
00457 dict["FEIMRC.a_tilt"] = feiexth.a_tilt;
00458 dict["FEIMRC.b_tilt"] = feiexth.b_tilt;
00459
00460 dict["FEIMRC.x_stage"] = feiexth.x_stage;
00461 dict["FEIMRC.y_stage"] = feiexth.y_stage;
00462 dict["FEIMRC.z_stage"] = feiexth.z_stage;
00463
00464 dict["FEIMRC.x_shift"] = feiexth.x_shift;
00465 dict["FEIMRC.y_shift"] = feiexth.y_shift;
00466
00467 dict["FEIMRC.defocus"] = feiexth.defocus;
00468 dict["FEIMRC.exp_time"] = feiexth.exp_time;
00469 dict["FEIMRC.mean_int"] = feiexth.mean_int;
00470 dict["FEIMRC.tilt_axis"] = feiexth.tilt_axis;
00471
00472 dict["FEIMRC.pixel_size"] = feiexth.pixel_size;
00473 dict["FEIMRC.magnification"] = feiexth.magnification;
00474 dict["FEIMRC.ht"] = feiexth.ht;
00475 dict["FEIMRC.binning"] = feiexth.binning;
00476 dict["FEIMRC.appliedDefocus"] = feiexth.appliedDefocus;
00477
00478
00479
00480 EXITFUNC;
00481 return 0;
00482 }
00483
00484 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00485 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00486 {
00487 ENTERFUNC;
00488
00489
00490 if(image_index == -1) {
00491 image_index = 0;
00492 }
00493 if(image_index != 0) {
00494 throw ImageWriteException(filename, "MRC file does not support stack.");
00495 }
00496 check_write_access(rw_mode, image_index, 1);
00497 if (area) {
00498 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00499 EXITFUNC;
00500 return 0;
00501 }
00502
00503 int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00504 int nx = dict["nx"];
00505 int ny = dict["ny"];
00506 int nz = dict["nz"];
00507 is_ri = dict["is_complex_ri"];
00508
00509 bool opposite_endian = false;
00510
00511 if (!is_new_file) {
00512 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00513 opposite_endian = true;
00514 }
00515 #if 0
00516 if (new_mode != mrch.mode) {
00517 LOGERR("cannot write to different mode file %s", filename.c_str());
00518 return 1;
00519 }
00520 #endif
00521 portable_fseek(mrcfile, 0, SEEK_SET);
00522 }
00523 else {
00524 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00525 mrch.mapc = 1;
00526 mrch.mapr = 2;
00527 mrch.maps = 3;
00528 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00529 }
00530
00531 if(nz<=1 && dict.has_key("xform.projection") && !dict.has_key("UCSF.chimera")) {
00532 Transform * t = dict["xform.projection"];
00533 Dict d = t->get_params("imagic");
00534 mrch.alpha = d["alpha"];
00535 mrch.beta = d["beta"];
00536 mrch.gamma = d["gamma"];
00537 mrch.xorigin = d["tx"];
00538 mrch.yorigin = d["ty"];
00539 mrch.zorigin = d["tz"];
00540 if(t) {delete t; t=0;}
00541 }
00542 else if(nz>1 && dict.has_key("xform.align3d") && !dict.has_key("UCSF.chimera")) {
00543 Transform * t = dict["xform.align3d"];
00544 Dict d = t->get_params("imagic");
00545 mrch.alpha = d["alpha"];
00546 mrch.beta = d["beta"];
00547 mrch.gamma = d["gamma"];
00548 mrch.xorigin = d["tx"];
00549 mrch.yorigin = d["ty"];
00550 mrch.zorigin = d["tz"];
00551 if(t) {delete t; t=0;}
00552 }
00553
00554 if(dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00555 mrch.xorigin = (float)dict["origin_x"];
00556 mrch.yorigin = (float)dict["origin_y"];
00557
00558 if (is_new_file) {
00559 mrch.zorigin = (float)dict["origin_z"];
00560 }
00561 else {
00562 mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00563 }
00564 }
00565
00566 if (dict.has_key("MRC.nlabels")) {
00567 mrch.nlabels = dict["MRC.nlabels"];
00568 }
00569
00570 for (int i = 0; i < MRC_NUM_LABELS; i++) {
00571 char label[32];
00572 #ifdef _WIN32
00573 _snprintf(label,31, "MRC.label%d", i);
00574 #else
00575 snprintf(label,31, "MRC.label%d", i);
00576 #endif //_WIN32
00577 if (dict.has_key(label)) {
00578 #ifdef _WIN32
00579 _snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00580 #else
00581 snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00582 #endif //_WIN32
00583 mrch.nlabels = i + 1;
00584 }
00585 }
00586
00587 if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00588 #ifdef _WIN32
00589 _snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00590 #else
00591 snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00592 #endif //_WIN32
00593 mrch.nlabels++;
00594 }
00595
00596 mrch.labels[mrch.nlabels][0] = '\0';
00597 mrch.mode = new_mode;
00598
00599 if (is_complex_mode()) {
00600 mrch.nx = nx / 2;
00601 }
00602 else {
00603 mrch.nx = nx;
00604 }
00605 mrch.ny = ny;
00606
00607 if (is_new_file) {
00608 mrch.nz = nz;
00609 }
00610 else if (image_index >= mrch.nz) {
00611 mrch.nz = image_index + 1;
00612 }
00613
00614 mrch.ispg = dict.has_key("MRC.ispg") ? (int)dict["MRC.ispg"] : 0;
00615 mrch.nsymbt = 0;
00616 mrch.amin = dict["minimum"];
00617 mrch.amax = dict["maximum"];
00618 mrch.amean = dict["mean"];
00619 mrch.rms = dict["sigma"];
00620
00623
00624
00625
00626
00627 mrch.mx = nx;
00628
00629
00630
00631
00632
00633 mrch.my = ny;
00634
00635
00636
00637
00638
00639 mrch.mz = nz;
00640
00641
00642 mrch.xlen = mrch.mx * (float) dict["apix_x"];
00643 mrch.ylen = mrch.my * (float) dict["apix_y"];
00644 mrch.zlen = mrch.mz * (float) dict["apix_z"];
00645
00646 if(dict.has_key("MRC.nxstart")) {
00647 mrch.nxstart = dict["MRC.nxstart"];
00648 }
00649 else {
00650 mrch.nxstart = -nx / 2;
00651 }
00652 if(dict.has_key("MRC.nystart")) {
00653 mrch.nystart = dict["MRC.nystart"];
00654 }
00655 else {
00656 mrch.nystart = -ny / 2;
00657 }
00658 if(dict.has_key("MRC.nzstart")) {
00659 mrch.nzstart = dict["MRC.nzstart"];
00660 }
00661 else {
00662 mrch.nzstart = -nz / 2;
00663 }
00664
00665 strncpy(mrch.map,"MAP ",4);
00666 mrch.machinestamp = generate_machine_stamp();
00667
00668 MrcHeader mrch2 = mrch;
00669
00670 if (opposite_endian || !use_host_endian) {
00671 swap_header(mrch2);
00672 }
00673
00674 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00675 throw ImageWriteException(filename, "MRC header");
00676 }
00677
00678 mode_size = get_mode_size(mrch.mode);
00679 is_new_file = false;
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689 EXITFUNC;
00690 return 0;
00691 }
00692
00693 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool )
00694 {
00695 ENTERFUNC;
00696
00697 if(!isFEI) {
00698
00699 image_index = 0;
00700 }
00701
00702 if(is_transpose && area!=0) {
00703 printf("Warning: This image dimension is in (y,x,z), region I/O not supported, return the whole image instead.");
00704 }
00705
00706 check_read_access(image_index, rdata);
00707
00708 if (area && is_complex_mode()) {
00709 LOGERR("Error: cannot read a region of a complex image.");
00710 return 1;
00711 }
00712
00713 unsigned char *cdata = (unsigned char *) rdata;
00714 short *sdata = (short *) rdata;
00715 unsigned short *usdata = (unsigned short *) rdata;
00716
00717 size_t size = 0;
00718 int xlen = 0, ylen = 0, zlen = 0;
00719 if(isFEI) {
00720 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
00721 portable_fseek(mrcfile, sizeof(MrcHeader)+feimrch.next, SEEK_SET);
00722
00723 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00724 image_index, mode_size,
00725 feimrch.nx, feimrch.ny, feimrch.nz, area);
00726
00727 EMUtil::get_region_dims(area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
00728
00729 size = (size_t)xlen * ylen * zlen;
00730 }
00731 else {
00732 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00733 portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00734
00735 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00736 image_index, mode_size,
00737 mrch.nx, mrch.ny, mrch.nz, area);
00738
00739 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00740
00741 size = (size_t)xlen * ylen * zlen;
00742 }
00743
00744 if (mrch.mode != MRC_UCHAR) {
00745 if (mode_size == sizeof(short)) {
00746 become_host_endian < short >(sdata, size);
00747 }
00748 else if (mode_size == sizeof(float)) {
00749 become_host_endian < float >(rdata, size);
00750 }
00751 }
00752
00753 if (mrch.mode == MRC_UCHAR) {
00754 for (size_t i = 0; i < size; ++i) {
00755 size_t j = size - 1 - i;
00756
00757 rdata[j] = static_cast < float >(cdata[j]);
00758 }
00759 }
00760 else if (mrch.mode == MRC_SHORT ) {
00761 for (size_t i = 0; i < size; ++i) {
00762 size_t j = size - 1 - i;
00763 rdata[j] = static_cast < float >(sdata[j]);
00764 }
00765 }
00766 else if (mrch.mode == MRC_USHORT) {
00767 for (size_t i = 0; i < size; ++i) {
00768 size_t j = size - 1 - i;
00769 rdata[j] = static_cast < float >(usdata[j]);
00770 }
00771 }
00772
00773 if(is_transpose) {
00774 transpose(rdata, xlen, ylen, zlen);
00775 }
00776
00777 if (is_complex_mode()) {
00778 if(!is_ri) Util::ap2ri(rdata, size);
00779 Util::flip_complex_phase(rdata, size);
00780 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00781 }
00782 EXITFUNC;
00783 return 0;
00784 }
00785
00786 int MrcIO::write_data(float *data, int image_index, const Region* area,
00787 EMUtil::EMDataType, bool use_host_endian)
00788 {
00789 ENTERFUNC;
00790
00791 image_index = 0;
00792 check_write_access(rw_mode, image_index, 1, data);
00793 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00794
00795 int nx, ny, nz;
00796 if(!area) {
00797 nx = mrch.nx;
00798 ny = mrch.ny;
00799 nz = mrch.nz;
00800 }
00801 else {
00802 nx = area->get_width();
00803 ny = area->get_height();
00804 nz = area->get_depth();
00805 }
00806 size_t size = (size_t)nx * ny * nz;
00807
00808 if (is_complex_mode()) {
00809 nx *= 2;
00810 if (!is_ri) {
00811 Util::ap2ri(data, size);
00812 is_ri = 1;
00813 }
00814 Util::flip_complex_phase(data, size);
00815 Util::rotate_phase_origin(data, nx, ny, nz);
00816 }
00817
00818 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00819
00820 if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian) {
00821 if (mrch.mode != MRC_UCHAR) {
00822 if (mode_size == sizeof(short)) {
00823 ByteOrder::swap_bytes((short*) data, size);
00824 }
00825 else if (mode_size == sizeof(float)) {
00826 ByteOrder::swap_bytes((float*) data, size);
00827 }
00828 }
00829 }
00830 mode_size = get_mode_size(mrch.mode);
00831
00832
00833
00834
00835 void * ptr_data = data;
00836
00837 float rendermin = 0.0f;
00838 float rendermax = 0.0f;
00839 EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, nz);
00840
00841 unsigned char *cdata = 0;
00842 short *sdata = 0;
00843 unsigned short *usdata = 0;
00844 if (mrch.mode == MRC_UCHAR) {
00845 cdata = new unsigned char[size];
00846 for (size_t i = 0; i < size; ++i) {
00847 if(data[i] <= rendermin) {
00848 cdata[i] = 0;
00849 }
00850 else if(data[i] >= rendermax){
00851 cdata[i] = UCHAR_MAX;
00852 }
00853 else {
00854 cdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
00855 }
00856 }
00857 ptr_data = cdata;
00858 update_stat((void *)cdata);
00859 }
00860 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00861 sdata = new short[size];
00862 for (size_t i = 0; i < size; ++i) {
00863 if(data[i] <= rendermin) {
00864 sdata[i] = SHRT_MIN;
00865 }
00866 else if(data[i] >= rendermax) {
00867 sdata[i] = SHRT_MAX;
00868 }
00869 else {
00870 sdata[i]=(short)(((data[i]-rendermin)/(rendermax-rendermin))*(SHRT_MAX-SHRT_MIN) - SHRT_MAX);
00871 }
00872 }
00873 ptr_data = sdata;
00874 update_stat((void *)sdata);
00875 }
00876 else if (mrch.mode == MRC_USHORT) {
00877 usdata = new unsigned short[size];
00878 for (size_t i = 0; i < size; ++i) {
00879 if(data[i] <= rendermin) {
00880 usdata[i] = 0;
00881 }
00882 else if(data[i] >= rendermax) {
00883 usdata[i] = USHRT_MAX;
00884 }
00885 else {
00886 usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
00887 }
00888 }
00889 ptr_data = usdata;
00890 update_stat((void *)usdata);
00891 }
00892
00893
00894
00895
00896 EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
00897 mode_size, mrch.nx, mrch.ny, mrch.nz, area);
00898
00899 if(cdata) {delete [] cdata; cdata=0;}
00900 if(sdata) {delete [] sdata; sdata=0;}
00901 if(usdata) {delete [] usdata; usdata=0;}
00902
00903 #if 0
00904 int row_size = nx * get_mode_size(mrch.mode);
00905 int sec_size = nx * ny;
00906
00907 unsigned char *cbuf = new unsigned char[row_size];
00908 unsigned short *sbuf = (unsigned short *) cbuf;
00909
00910 for (int i = 0; i < nz; i++) {
00911 int i2 = i * sec_size;
00912 for (int j = 0; j < ny; j++) {
00913 int k = i2 + j * nx;
00914 void *pbuf = 0;
00915
00916 switch (mrch.mode) {
00917 case MRC_UCHAR:
00918 for (int l = 0; l < nx; l++) {
00919 cbuf[l] = static_cast < unsigned char >(data[k + l]);
00920 }
00921 pbuf = cbuf;
00922 fwrite(cbuf, row_size, 1, mrcfile);
00923 break;
00924
00925 case MRC_SHORT:
00926 case MRC_SHORT_COMPLEX:
00927 for (int l = 0; l < nx; l++) {
00928 sbuf[l] = static_cast < short >(data[k + l]);
00929 }
00930 pbuf = sbuf;
00931 fwrite(sbuf, row_size, 1, mrcfile);
00932 break;
00933
00934 case MRC_USHORT:
00935 for (int l = 0; l < nx; l++) {
00936 sbuf[l] = static_cast < unsigned short >(data[k + l]);
00937 }
00938 pbuf = sbuf;
00939 fwrite(sbuf, row_size, 1, mrcfile);
00940 break;
00941
00942 case MRC_FLOAT:
00943 case MRC_FLOAT_COMPLEX:
00944 pbuf = &data[k];
00945 break;
00946 }
00947 if (pbuf) {
00948 fwrite(pbuf, row_size, 1, mrcfile);
00949 }
00950 }
00951 }
00952
00953 if(cbuf)
00954 {
00955 delete[]cbuf;
00956 cbuf = 0;
00957 }
00958 #endif
00959
00960 EXITFUNC;
00961 return 0;
00962 }
00963
00964 void MrcIO::update_stat(void* data)
00965 {
00966 size_t size = mrch.nx * mrch.ny * mrch.nz;
00967 float v = 0.0f;
00968 double sum = 0.0;
00969 double square_sum = 0.0;
00970 double mean = 0.0;
00971 float min, max;
00972
00973 unsigned char * cdata = 0;
00974 short * sdata = 0;
00975 unsigned short * usdata = 0;
00976
00977 if (mrch.mode == MRC_UCHAR) {
00978 max = 0.0f;
00979 min = UCHAR_MAX;
00980 cdata = (unsigned char *)data;
00981
00982 for (size_t i = 0; i < size; ++i) {
00983 v = (float)(cdata[i]);
00984 #ifdef _WIN32
00985 max = _cpp_max(max,v);
00986 min = _cpp_min(min,v);
00987 #else
00988 max=std::max<float>(max,v);
00989 min=std::min<float>(min,v);
00990 #endif //_WIN32
00991
00992 sum += v;
00993 square_sum += v * v;
00994 }
00995 }
00996 else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00997 max = (float)SHRT_MIN;
00998 min = (float)SHRT_MAX;
00999 sdata = (short *)data;
01000
01001 for (size_t i = 0; i < size; ++i) {
01002 v = (float)(sdata[i]);
01003 #ifdef _WIN32
01004 max = _cpp_max(max,v);
01005 min = _cpp_min(min,v);
01006 #else
01007 max=std::max<float>(max,v);
01008 min=std::min<float>(min,v);
01009 #endif //_WIN32
01010
01011 sum += v;
01012 square_sum += v * v;
01013 }
01014 }
01015 else if (mrch.mode == MRC_USHORT) {
01016 max = 0.0f;
01017 min = (float)USHRT_MAX;
01018 usdata = (unsigned short*)data;
01019
01020 for (size_t i = 0; i < size; ++i) {
01021 v = (float)(usdata[i]);
01022 #ifdef _WIN32
01023 max = _cpp_max(max,v);
01024 min = _cpp_min(min,v);
01025 #else
01026 max=std::max<float>(max,v);
01027 min=std::min<float>(min,v);
01028 #endif //_WIN32
01029
01030 sum += v;
01031 square_sum += v * v;
01032 }
01033 }
01034 else {
01035 throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
01036 }
01037
01038 mean = sum/size;
01039 #ifdef _WIN32
01040 float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / size)/(size-1)));
01041 #else
01042 float sigma = (float)std::sqrt(std::max<float>(0.0,(square_sum - sum*sum / size)/(size-1)));
01043 #endif //_WIN32
01044
01045
01046 mrch.amin = min;
01047 mrch.amax = max;
01048 mrch.amean = (float)mean;
01049 mrch.rms = sigma;
01050
01051 MrcHeader mrch2 = mrch;
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065
01066
01067
01068 portable_fseek(mrcfile, 0, SEEK_SET);
01069
01070 if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
01071 throw ImageWriteException(filename, "MRC header");
01072 }
01073
01074 portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
01075 }
01076
01077 bool MrcIO::is_complex_mode()
01078 {
01079 init();
01080 if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
01081 return true;
01082 }
01083 return false;
01084 }
01085
01086
01087 int MrcIO::read_ctf(Ctf & ctf, int)
01088 {
01089 ENTERFUNC;
01090 init();
01091 size_t n = strlen(CTF_MAGIC);
01092
01093 int err = 1;
01094 if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
01095 err = ctf.from_string(string(&mrch.labels[0][n]));
01096 }
01097 EXITFUNC;
01098 return err;
01099 }
01100
01101 void MrcIO::write_ctf(const Ctf & ctf, int)
01102 {
01103 ENTERFUNC;
01104 init();
01105
01106 string ctf_str = ctf.to_string();
01107 #ifdef _WIN32
01108 _snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01109 #else
01110 snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01111 #endif //_WIN32
01112 rewind(mrcfile);
01113
01114 if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
01115 throw ImageWriteException(filename, "write CTF info to header failed");
01116 }
01117 EXITFUNC;
01118 }
01119
01120 void MrcIO::flush()
01121 {
01122 fflush(mrcfile);
01123 }
01124
01125
01126 int MrcIO::get_mode_size(int mm)
01127 {
01128 MrcIO::MrcMode m = static_cast < MrcMode > (mm);
01129
01130 int msize = 0;
01131 switch (m) {
01132 case MRC_UCHAR:
01133 msize = sizeof(char);
01134 break;
01135 case MRC_SHORT:
01136 case MRC_USHORT:
01137 case MRC_SHORT_COMPLEX:
01138 msize = sizeof(short);
01139 break;
01140 case MRC_FLOAT:
01141 case MRC_FLOAT_COMPLEX:
01142 msize = sizeof(float);
01143 break;
01144 default:
01145 msize = 0;
01146 }
01147
01148 return msize;
01149 }
01150
01151 int MrcIO::to_em_datatype(int m)
01152 {
01153 EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
01154
01155 switch (m) {
01156 case MRC_UCHAR:
01157 e = EMUtil::EM_UCHAR;
01158 break;
01159 case MRC_SHORT:
01160 e = EMUtil::EM_SHORT;
01161 break;
01162 case MRC_USHORT:
01163 e = EMUtil::EM_USHORT;
01164 break;
01165 case MRC_SHORT_COMPLEX:
01166 e = EMUtil::EM_SHORT_COMPLEX;
01167 break;
01168 case MRC_FLOAT:
01169 e = EMUtil::EM_FLOAT;
01170 break;
01171 case MRC_FLOAT_COMPLEX:
01172 e = EMUtil::EM_FLOAT_COMPLEX;
01173 break;
01174 default:
01175 e = EMUtil::EM_UNKNOWN;
01176 }
01177 return e;
01178 }
01179
01180
01181 int MrcIO::to_mrcmode(int e, int is_complex)
01182 {
01183 MrcMode m = MRC_UNKNOWN;
01184 EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
01185
01186 switch (em_type) {
01187 case EMUtil::EM_UCHAR:
01188 m = MRC_UCHAR;
01189 break;
01190 case EMUtil::EM_USHORT:
01191 if (is_complex) {
01192 m = MRC_SHORT_COMPLEX;
01193 }
01194 else {
01195 m = MRC_USHORT;
01196 }
01197 break;
01198 case EMUtil::EM_SHORT:
01199 if (is_complex) {
01200 m = MRC_SHORT_COMPLEX;
01201 }
01202 else {
01203 m = MRC_SHORT;
01204 }
01205 break;
01206 case EMUtil::EM_SHORT_COMPLEX:
01207 case EMUtil::EM_USHORT_COMPLEX:
01208 m = MRC_SHORT_COMPLEX;
01209 break;
01210 case EMUtil::EM_CHAR:
01211 case EMUtil::EM_INT:
01212 case EMUtil::EM_UINT:
01213 case EMUtil::EM_FLOAT:
01214 if (is_complex) {
01215 m = MRC_FLOAT_COMPLEX;
01216 }
01217 else {
01218 m = MRC_FLOAT;
01219 }
01220 break;
01221 case EMUtil::EM_FLOAT_COMPLEX:
01222 m = MRC_FLOAT_COMPLEX;
01223 break;
01224 default:
01225 m = MRC_FLOAT;
01226 }
01227
01228 return m;
01229 }
01230
01231
01232
01233 int MrcIO::generate_machine_stamp()
01234 {
01235 int stamp = 0;
01236 char *p = (char *) (&stamp);
01237
01238 if (ByteOrder::is_host_big_endian()) {
01239 p[0] = 0x11;
01240 p[1] = 0x11;
01241 p[2] = 0;
01242 p[3] = 0;
01243 }
01244 else {
01245 p[0] = 0x44;
01246 p[1] = 0x41;
01247 p[2] = 0;
01248 p[3] = 0;
01249 }
01250 return stamp;
01251 }
01252
01253 void MrcIO::swap_header(MrcHeader& mrch)
01254 {
01255 ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
01256 ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
01257 }
01258
01259 int MrcIO::get_nimg()
01260 {
01261 init();
01262
01263 return 1;
01264 }
01265
01266 int MrcIO::transpose(float *data, int xlen, int ylen, int zlen) const
01267 {
01268 float * tmp = new float[xlen*ylen];
01269
01270 for(size_t z=0; z<(size_t)zlen; ++z) {
01271 for(size_t y=0; y<(size_t)ylen; ++y) {
01272 for(size_t x=0; x<(size_t)xlen; ++x) {
01273 tmp[x*ylen+y] = data[z*xlen*ylen+y*xlen+x];
01274 }
01275 }
01276 std::copy(tmp, tmp+xlen*ylen, data+z*xlen*ylen);
01277 }
01278
01279 delete [] tmp;
01280
01281 return 0;
01282 }