EMAN2
mrcio.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include <cstring>
00037 #include <climits>
00038 
00039 #include <sys/stat.h>
00040 
00041 #include "mrcio.h"
00042 #include "portable_fileio.h"
00043 #include "geometry.h"
00044 #include "util.h"
00045 #include "ctf.h"
00046 #include "transform.h"
00047 
00048 using namespace EMAN;
00049 
00050 const char *MrcIO::CTF_MAGIC = "!-";
00051 const char *MrcIO::SHORT_CTF_MAGIC = "!$";
00052 
00053 MrcIO::MrcIO(const string & mrc_filename, IOMode rw)
00054 :       filename(mrc_filename), rw_mode(rw), mrcfile(0), mode_size(0),
00055                 isFEI(false), is_ri(0), is_new_file(false), initialized(false),
00056                 is_transpose(false), is_stack(false), stack_size(1)
00057 {
00058         memset(&mrch, 0, sizeof(MrcHeader));
00059         is_big_endian = ByteOrder::is_host_big_endian();
00060 }
00061 
00062 MrcIO::~MrcIO()
00063 {
00064         if (mrcfile) {
00065                 fclose(mrcfile);
00066                 mrcfile = NULL;
00067         }
00068 }
00069 
00070 void MrcIO::init()
00071 {
00072         ENTERFUNC;
00073 
00074         if (initialized) {
00075                 return;
00076         }
00077 
00078         setbuf (stdout, NULL);
00079 
00080         IOMode rwmode;
00081 
00082         if (rw_mode == WRITE_ONLY) {
00083                 rwmode = READ_WRITE;
00084         }
00085         else {
00086                 rwmode = rw_mode;
00087         }
00088 
00089         int error_type;
00090         struct stat status;
00091 
00092         error_type = stat(filename.c_str(), & status);
00093         is_new_file = (error_type != 0);
00094 
00095         initialized = true;
00096 
00097 //      mrcfile = sfopen(filename, rwmode, &is_new_file);
00098         mrcfile = sfopen(filename, rwmode, NULL);
00099 
00100         string ext = Util::get_filename_ext(filename);
00101 
00102         if (ext != "") {
00103                 if (ext == "raw"   ||  ext == "RAW") {
00104                         isFEI = true;
00105                 }
00106 
00107                 if (ext == "mrcs"  ||  ext == "MRCS") {
00108                         is_stack = true;
00109                 }
00110         }
00111 
00112         if (! is_new_file) {
00113                 if (fread(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00114                         throw ImageReadException(filename, "MRC header");
00115                 }
00116 
00117                 if (! is_valid(&mrch)) {
00118                         throw ImageReadException(filename, "invalid MRC");
00119                 }
00120 
00121                 is_big_endian = ByteOrder::is_data_big_endian(&mrch.nz);
00122 
00123                 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00124                         swap_header(mrch);
00125                 }
00126 
00127                 // become_host_endian((int *) &mrch, NUM_4BYTES_PRE_MAP);
00128                 // become_host_endian((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
00129 
00130                 mode_size = get_mode_size(mrch.mode);
00131 
00132                 if (is_complex_mode()) {
00133                         is_ri = 1;
00134                 }
00135 
00136                 if (mrch.nxstart != 0 || mrch.nystart != 0 || mrch.nzstart != 0) {
00137                         LOGWARN("nx/ny/nz start not zero");
00138                 }
00139 
00140                 if (is_complex_mode()) {
00141                         mrch.nx *= 2;
00142                 }
00143 
00144                 if (mrch.xlen == 0) {
00145                         mrch.xlen = 1.0;
00146                 }
00147 
00148                 if (mrch.ylen == 0) {
00149                         mrch.ylen = 1.0;
00150                 }
00151 
00152                 if (mrch.zlen == 0) {
00153                         mrch.zlen = 1.0;
00154                 }
00155 
00156                 if (mrch.nlabels > 0) {
00157                         if (string(mrch.labels[0],3) == "Fei") {
00158                                 isFEI = true;
00159                         }
00160                 }
00161 
00162                 if (mrch.mapc == 2 && mrch.mapr == 1) {
00163                         is_transpose = true;
00164                 }
00165 
00166                 if (is_stack) {
00167                         stack_size = mrch.nz;
00168                         mrch.nz = 1;
00169                 }
00170         }
00171 
00172         EXITFUNC;
00173 }
00174 
00175 bool MrcIO::is_image_big_endian()
00176 {
00177         init();
00178 
00179         return is_big_endian;
00180 }
00181 
00182 bool MrcIO::is_valid(const void *first_block, off_t file_size)
00183 {
00184         ENTERFUNC;
00185 
00186         if (!first_block) {
00187                 return false;
00188         }
00189 
00190         const int *data = static_cast < const int *>(first_block);
00191         int nx = data[0];
00192         int ny = data[1];
00193         int nz = data[2];
00194         int mrcmode = data[3];
00195         int nsymbt = data[23];  // this field specify the extra bytes for symmetry information
00196 
00197         bool data_big_endian = ByteOrder::is_data_big_endian(&nz);
00198 
00199         if (data_big_endian != ByteOrder::is_host_big_endian()) {
00200                 ByteOrder::swap_bytes(&nx);
00201                 ByteOrder::swap_bytes(&ny);
00202                 ByteOrder::swap_bytes(&nz);
00203                 ByteOrder::swap_bytes(&mrcmode);
00204                 ByteOrder::swap_bytes(&nsymbt);
00205         }
00206 
00207         if (mrcmode == MRC_SHORT_COMPLEX || mrcmode == MRC_FLOAT_COMPLEX) {
00208                 nx *= 2;
00209         }
00210 
00211         const int max_dim = 1 << 20;
00212 
00213         if ((mrcmode >= MRC_UCHAR && mrcmode < MRC_UNKNOWN) &&
00214                 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim)
00215                  && (nz > 0 && nz < max_dim)) {
00216 
00217 //#ifndef SPIDERMRC // Spider MRC files don't satisfy the following test
00218 
00219                 if (file_size > 0) {
00220                         off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz *
00221                                 (off_t)get_mode_size(mrcmode) +
00222                                 (off_t)sizeof(MrcHeader) + nsymbt;
00223 
00224                         if (file_size == file_size1) {
00225                                 return true;
00226                         }
00227 
00228 //                      return false;
00229 
00230          // when size doesn't match, print error message instead of make it fail
00231 
00232                         LOGWARN("image size check fails, still try to read it...");
00233                 }
00234                 else {
00235                         return true;
00236                 }
00237 
00238 //#endif // SPIDERMRC
00239 
00240                 return true;
00241         }
00242 
00243         EXITFUNC;
00244 
00245         return false;
00246 }
00247 
00248 int MrcIO::read_header(Dict & dict, int image_index, const Region * area, bool is_3d)
00249 {
00250         init();
00251 
00252         if (isFEI) {
00253                 return read_fei_header(dict, image_index, area, is_3d);
00254         }
00255         else {
00256                 return read_mrc_header(dict, image_index, area, is_3d);
00257         }
00258 }
00259 
00260 int MrcIO::read_mrc_header(Dict & dict, int image_index, const Region * area, bool)
00261 {
00262         ENTERFUNC;
00263 
00264         if (image_index < 0) {
00265                 image_index = 0;
00266         }
00267 
00268         if (image_index != 0  &&  ! is_stack) {
00269                 throw ImageReadException(filename,
00270                         "no stack allowed for MRC image. For take 2D slice out of 3D image, "
00271                         "read the 3D image first, then use get_clip().");
00272         }
00273 
00274         check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00275 
00276         int xlen = 0, ylen = 0, zlen = 0;
00277 
00278         EMUtil::get_region_dims(area, mrch.nx, & xlen, mrch.ny, & ylen, mrch.nz, & zlen);
00279 
00280         dict["nx"] = xlen;
00281         dict["ny"] = ylen;
00282         dict["nz"] = zlen;
00283         dict["MRC.nx"] = mrch.nx;
00284         dict["MRC.ny"] = mrch.ny;
00285         dict["MRC.nz"] = mrch.nz;
00286 
00287         dict["datatype"] = to_em_datatype(mrch.mode);
00288 
00289         dict["MRC.nxstart"] = mrch.nxstart;
00290         dict["MRC.nystart"] = mrch.nystart;
00291         dict["MRC.nzstart"] = mrch.nzstart;
00292 
00293         dict["MRC.mx"] = mrch.mx;
00294         dict["MRC.my"] = mrch.my;
00295         dict["MRC.mz"] = mrch.mz;
00296 
00297         dict["MRC.xlen"] = mrch.xlen;
00298         dict["MRC.ylen"] = mrch.ylen;
00299         dict["MRC.zlen"] = mrch.zlen;
00300 
00301         dict["MRC.alpha"] = mrch.alpha;
00302         dict["MRC.beta"] = mrch.beta;
00303         dict["MRC.gamma"] = mrch.gamma;
00304 
00305         dict["MRC.mapc"] = mrch.mapc;
00306         dict["MRC.mapr"] = mrch.mapr;
00307         dict["MRC.maps"] = mrch.maps;
00308 
00309         dict["MRC.minimum"] = mrch.amin;
00310         dict["MRC.maximum"] = mrch.amax;
00311         dict["MRC.mean"] = mrch.amean;
00312         dict["minimum"] = mrch.amin;
00313         dict["maximum"] = mrch.amax;
00314         dict["mean"] = mrch.amean;
00315 
00316         dict["MRC.ispg"] = mrch.ispg;
00317         dict["MRC.nsymbt"] = mrch.nsymbt;
00318 
00319         float apx = mrch.xlen / mrch.mx;
00320         float apy = mrch.ylen / mrch.my;
00321         float apz = mrch.zlen / mrch.mz;
00322 
00323         if (apx > 1000 || apx < 0.01) {
00324                 dict["apix_x"] = 1.0f;
00325         }
00326         else {
00327                 dict["apix_x"] = apx;
00328         }
00329 
00330         if (apy > 1000 || apy < 0.01) {
00331                 dict["apix_y"] = 1.0f;
00332         }
00333         else {
00334                 dict["apix_y"] = apy;
00335         }
00336 
00337         if (apz > 1000 || apz < 0.01) {
00338                 dict["apix_z"] = 1.0f;
00339         }
00340         else {
00341                 dict["apix_z"] = apz;
00342         }
00343 
00344         if (area) {
00345                 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
00346                 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
00347 
00348                 if (area->get_ndim() == 3 && mrch.nz > 1) {
00349                         dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
00350                 }
00351                 else {
00352                         dict["origin_z"] = mrch.zorigin;
00353                 }
00354         }
00355         else {
00356                 dict["origin_x"] = mrch.xorigin;
00357                 dict["origin_y"] = mrch.yorigin;
00358                 dict["origin_z"] = mrch.zorigin;
00359         }
00360 
00361         if (is_complex_mode()) {
00362                 dict["is_complex"] = 1;
00363                 dict["is_complex_ri"] = 1;
00364         }
00365 
00366         dict["MRC.machinestamp"] = mrch.machinestamp;
00367 
00368         dict["MRC.rms"] = mrch.rms;
00369         dict["sigma"] = mrch.rms;
00370         dict["MRC.nlabels"] = mrch.nlabels;
00371 
00372         for (int i = 0; i < mrch.nlabels; i++) {
00373                 char label[32];
00374                 sprintf(label, "MRC.label%d", i);
00375                 dict[string(label)] = string(mrch.labels[i],80);
00376         }
00377 
00378         EMAN1Ctf ctf_;
00379 
00380         if (read_ctf(ctf_) == 0) {
00381                 vector<float> vctf = ctf_.to_vector();
00382                 dict["ctf"] = vctf;
00383         }
00384 
00385         if (is_transpose) {
00386                 dict["nx"] = ylen;
00387                 dict["ny"] = xlen;
00388                 dict["MRC.nx"] = mrch.ny;
00389                 dict["MRC.ny"] = mrch.nx;
00390                 dict["MRC.mx"] = mrch.my;
00391                 dict["MRC.my"] = mrch.mx;
00392                 dict["apix_x"] = mrch.ylen / mrch.my;
00393                 dict["apix_y"] = mrch.xlen / mrch.mx;
00394                 dict["origin_x"] = mrch.yorigin;
00395                 dict["origin_y"] = mrch.xorigin;
00396                 dict["MRC.nxstart"] = mrch.nystart;
00397                 dict["MRC.nystart"] = mrch.nxstart;
00398         }
00399 
00400         Transform * trans = new Transform();
00401 
00402         if (is_transpose) {
00403                 trans->set_trans(mrch.nystart, mrch.nxstart, mrch.nzstart);
00404                 trans->set_rotation(Dict("type", "imagic", "alpha", mrch.alpha,
00405                                                                                  "beta", mrch.beta, "gamma", mrch.gamma));
00406         }
00407         else {
00408                 trans->set_trans(mrch.nxstart, mrch.nystart, mrch.nzstart);
00409                 trans->set_rotation(Dict("type", "imagic", "alpha", mrch.alpha,
00410                                                                                  "beta", mrch.beta, "gamma", mrch.gamma));
00411         }
00412         
00413         if (zlen <= 1) {
00414                 dict["xform.projection"] = trans;
00415         }
00416         else {
00417                 dict["xform.align3d"] = trans;
00418         }
00419 
00420         if (trans) {delete trans; trans = NULL;}
00421 
00422         EXITFUNC;
00423 
00424         return 0;
00425 }
00426 
00427 int MrcIO::read_fei_header(Dict & dict, int image_index, const Region * area, bool)
00428 {
00429         ENTERFUNC;
00430 
00431         if(image_index < 0) {
00432                 image_index = 0;
00433         }
00434 
00435         init();
00436 
00437         check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
00438 
00439         int xlen = 0, ylen = 0, zlen = 0;
00440 
00441         EMUtil::get_region_dims(area, feimrch.nx, & xlen, feimrch.ny, & ylen, feimrch.nz, & zlen);
00442 
00443         dict["nx"] = xlen;
00444         dict["ny"] = ylen;
00445         dict["nz"] = zlen;
00446         dict["FEIMRC.nx"] = feimrch.nx;
00447         dict["FEIMRC.ny"] = feimrch.ny;
00448         dict["FEIMRC.nz"] = feimrch.nz;
00449 
00450         dict["datatype"] = to_em_datatype(feimrch.mode);        //=1, FEI-MRC file always use short for data type
00451 
00452         dict["FEIMRC.nxstart"] = feimrch.nxstart;
00453         dict["FEIMRC.nystart"] = feimrch.nystart;
00454         dict["FEIMRC.nzstart"] = feimrch.nzstart;
00455 
00456         dict["FEIMRC.mx"] = feimrch.mx;
00457         dict["FEIMRC.my"] = feimrch.my;
00458         dict["FEIMRC.mz"] = feimrch.mz;
00459 
00460         dict["FEIMRC.xlen"] = feimrch.xlen;
00461         dict["FEIMRC.ylen"] = feimrch.ylen;
00462         dict["FEIMRC.zlen"] = feimrch.zlen;
00463 
00464         dict["FEIMRC.alpha"] = feimrch.alpha;
00465         dict["FEIMRC.beta"] = feimrch.beta;
00466         dict["FEIMRC.gamma"] = feimrch.gamma;
00467 
00468         dict["FEIMRC.mapc"] = feimrch.mapc;
00469         dict["FEIMRC.mapr"] = feimrch.mapr;
00470         dict["FEIMRC.maps"] = feimrch.maps;
00471 
00472         dict["FEIMRC.minimum"] = feimrch.amin;
00473         dict["FEIMRC.maximum"] = feimrch.amax;
00474         dict["FEIMRC.mean"] = feimrch.amean;
00475         dict["mean"] = feimrch.amean;
00476 
00477         dict["FEIMRC.ispg"] = feimrch.ispg;
00478         dict["FEIMRC.nsymbt"] = feimrch.nsymbt;
00479 
00480         dict["apix_x"] = feimrch.xlen / feimrch.mx;
00481         dict["apix_y"] = feimrch.ylen / feimrch.my;
00482         dict["apix_z"] = feimrch.zlen / feimrch.mz;
00483 
00484         dict["FEIMRC.next"] = feimrch.next;     //offset from end of header to the first dataset
00485         dict["FEIMRC.dvid"] = feimrch.dvid;
00486         dict["FEIMRC.numintegers"] = feimrch.numintegers;
00487         dict["FEIMRC.numfloats"] = feimrch.numfloats;
00488         dict["FEIMRC.sub"] = feimrch.sub;
00489         dict["FEIMRC.zfac"] = feimrch.zfac;
00490 
00491         dict["FEIMRC.min2"] = feimrch.min2;
00492         dict["FEIMRC.max2"] = feimrch.max2;
00493         dict["FEIMRC.min3"] = feimrch.min3;
00494         dict["FEIMRC.max3"] = feimrch.max3;
00495         dict["FEIMRC.min4"] = feimrch.min4;
00496         dict["FEIMRC.max4"] = feimrch.max4;
00497 
00498         dict["FEIMRC.idtype"] = feimrch.idtype;
00499         dict["FEIMRC.lens"] = feimrch.lens;
00500         dict["FEIMRC.nd1"] = feimrch.nd1;
00501         dict["FEIMRC.nd2"] = feimrch.nd2;
00502         dict["FEIMRC.vd1"] = feimrch.vd1;
00503         dict["FEIMRC.vd2"] = feimrch.vd2;
00504 
00505         for(int i=0; i<9; i++) {        // 9 tilt angles
00506                 char label[32];
00507                 sprintf(label, "MRC.tiltangles%d", i);
00508                 dict[string(label)] = feimrch.tiltangles[i];
00509         }
00510 
00511         dict["FEIMRC.zorg"] = feimrch.zorg;
00512         dict["FEIMRC.xorg"] = feimrch.xorg;
00513         dict["FEIMRC.yorg"] = feimrch.yorg;
00514 
00515         dict["FEIMRC.nlabl"] = feimrch.nlabl;
00516 
00517         for (int i = 0; i < feimrch.nlabl; i++) {
00518                 char label[32];
00519                 sprintf(label, "MRC.label%d", i);
00520                 dict[string(label)] = string(feimrch.labl[i], 80);
00521         }
00522 
00523         /* Read extended image header by specified image index */
00524 
00525         FeiMrcExtHeader feiexth;
00526 
00527         portable_fseek(mrcfile, sizeof(FeiMrcHeader)+sizeof(FeiMrcExtHeader)*image_index, SEEK_SET);
00528 
00529         if (fread(&feiexth, sizeof(FeiMrcExtHeader), 1, mrcfile) != 1) {
00530                 throw ImageReadException(filename, "FEI MRC extended header");
00531         }
00532 
00533         dict["FEIMRC.a_tilt"] = feiexth.a_tilt;
00534         dict["FEIMRC.b_tilt"] = feiexth.b_tilt;
00535 
00536         dict["FEIMRC.x_stage"] = feiexth.x_stage;
00537         dict["FEIMRC.y_stage"] = feiexth.y_stage;
00538         dict["FEIMRC.z_stage"] = feiexth.z_stage;
00539 
00540         dict["FEIMRC.x_shift"] = feiexth.x_shift;
00541         dict["FEIMRC.y_shift"] = feiexth.y_shift;
00542 
00543         dict["FEIMRC.defocus"] = feiexth.defocus;
00544         dict["FEIMRC.exp_time"] = feiexth.exp_time;
00545         dict["FEIMRC.mean_int"] = feiexth.mean_int;
00546         dict["FEIMRC.tilt_axis"] = feiexth.tilt_axis;
00547 
00548         dict["FEIMRC.pixel_size"] = feiexth.pixel_size;
00549         dict["FEIMRC.magnification"] = feiexth.magnification;
00550         dict["FEIMRC.ht"] = feiexth.ht;
00551         dict["FEIMRC.binning"] = feiexth.binning;
00552         dict["FEIMRC.appliedDefocus"] = feiexth.appliedDefocus;
00553 
00554         // remainder 16 4-byte floats not used
00555 
00556         EXITFUNC;
00557 
00558         return 0;
00559 }
00560 
00561 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00562                                                 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00563 {
00564         ENTERFUNC;
00565 
00566         string ext = Util::get_filename_ext(filename);
00567         is_stack = (ext == "mrcs"  ||  ext == "MRCS");
00568 
00569         bool append = (image_index == -1);
00570 
00571         if (image_index == -1) {
00572                 image_index = 0;
00573         }
00574 
00575         if (image_index != 0  &&  ! is_stack) {
00576                 throw ImageWriteException(filename, "MRC file does not support stack.");
00577         }
00578 
00579         int max_images = 0;
00580 
00581         if (! is_stack) {
00582                 max_images = 1;
00583         }
00584 
00585         check_write_access(rw_mode, image_index, max_images);
00586 
00587         if (area) {
00588                 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00589 
00590                 EXITFUNC;
00591 
00592                 return 0;
00593         }
00594 
00595         int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00596         int nx = dict["nx"];
00597         int ny = dict["ny"];
00598         int nz = dict["nz"];
00599 
00600         bool got_one_image = (nz > 1);
00601 
00602         is_ri =  dict["is_complex_ri"];
00603 
00604         bool opposite_endian = false;
00605 
00606         if (! is_new_file) {
00607                 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00608                         opposite_endian = true;
00609                 }
00610 
00611 #if 0
00612                 if (new_mode != mrch.mode) {
00613                         LOGERR("cannot write to different mode file %s", filename.c_str());
00614 
00615                         return 1;
00616                 }
00617 #endif
00618 
00619                 portable_fseek(mrcfile, 0, SEEK_SET);
00620         }
00621         else {
00622                 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00623                 mrch.mapc = 1;
00624                 mrch.mapr = 2;
00625                 mrch.maps = 3;
00626                 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00627         }
00628 
00629         if (nz <= 1 && dict.has_key("xform.projection") && ! dict.has_key("UCSF.chimera")) {
00630                 Transform * t = dict["xform.projection"];
00631                 Dict d = t->get_params("imagic");
00632                 mrch.alpha   = d["alpha"];
00633                 mrch.beta    = d["beta"];
00634                 mrch.gamma   = d["gamma"];
00635                 mrch.xorigin = d["tx"];
00636                 mrch.yorigin = d["ty"];
00637                 mrch.zorigin = d["tz"];
00638 
00639                 if (t) {delete t; t = NULL;}
00640         }
00641         else if (nz > 1 && dict.has_key("xform.align3d") && ! dict.has_key("UCSF.chimera")) {
00642                 Transform * t = dict["xform.align3d"];
00643                 Dict d = t->get_params("imagic");
00644                 mrch.alpha   = d["alpha"];
00645                 mrch.beta    = d["beta"];
00646                 mrch.gamma   = d["gamma"];
00647                 mrch.xorigin = d["tx"];
00648                 mrch.yorigin = d["ty"];
00649                 mrch.zorigin = d["tz"];
00650 
00651                 if (t) {delete t; t = NULL;}
00652         }
00653 
00654         if (dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00655                 mrch.xorigin = (float)dict["origin_x"];
00656                 mrch.yorigin = (float)dict["origin_y"];
00657 
00658                 if (is_new_file) {
00659                         mrch.zorigin = (float)dict["origin_z"];
00660                 }
00661                 else {
00662                         mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00663                 }
00664         }
00665 
00666         if (dict.has_key("MRC.nlabels")) {
00667                 mrch.nlabels = dict["MRC.nlabels"];
00668         }
00669 
00670         for (int i = 0; i < MRC_NUM_LABELS; i++) {
00671                 char label[32];
00672 
00673 #ifdef _WIN32
00674                 _snprintf(label,31, "MRC.label%d", i);
00675 #else
00676                  snprintf(label,31, "MRC.label%d", i);
00677 #endif  // _WIN32
00678 
00679                 if (dict.has_key(label)) {
00680 #ifdef _WIN32
00681                         _snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00682 #else
00683                          snprintf(&mrch.labels[i][0],80, "%s", (const char *) dict[label]);
00684 #endif  // _WIN32
00685 
00686                         mrch.nlabels = i + 1;
00687                 }
00688         }
00689 
00690         if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00691 #ifdef _WIN32
00692                 _snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00693 #else
00694                  snprintf(&mrch.labels[mrch.nlabels][0],79, "EMAN %s", Util::get_time_label().c_str());
00695 #endif  // _WIN32
00696                 mrch.nlabels++;
00697         }
00698 
00699         mrch.labels[mrch.nlabels][0] = '\0';
00700         mrch.mode = new_mode;
00701 
00702         if (is_complex_mode()) {
00703                 mrch.nx = nx / 2;
00704         }
00705         else {
00706                 mrch.nx = nx;
00707         }
00708 
00709         mrch.ny = ny;
00710 
00711         if (is_stack) {
00712                 if (got_one_image) {
00713                         stack_size = nz;
00714                         image_index = 0;
00715                 }
00716                 else if (is_new_file) {
00717                         stack_size = 1;
00718                         image_index = stack_size - 1;
00719                 }
00720                 else if (append) {
00721                         stack_size++;
00722                         image_index = stack_size - 1;
00723                 }
00724                 else if (image_index >= stack_size) {
00725                         stack_size = image_index + 1;
00726                 }
00727 
00728                 nz = stack_size;
00729                 mrch.nz = nz;
00730         }
00731         else {
00732                 if (is_new_file) {
00733                         mrch.nz = nz;
00734                 }
00735                 else if (append) {
00736                         nz++;
00737                         mrch.nz = nz;
00738                 }
00739                 else if (image_index >= nz) {
00740                         nz = image_index + 1;
00741                         mrch.nz = nz;
00742                 }
00743                 else {
00744                         mrch.nz = nz;
00745                 }
00746 
00747                 stack_size = 1;
00748         }
00749 
00750         mrch.ispg = dict.has_key("MRC.ispg") ? (int)dict["MRC.ispg"] : 0;
00751         mrch.nsymbt = 0;
00752         mrch.amin = dict["minimum"];
00753         mrch.amax = dict["maximum"];
00754         mrch.amean = dict["mean"];
00755         mrch.rms = dict["sigma"];
00756 
00760 //      if(dict.has_key("MRC.mx")) {
00761 //              mrch.mx = dict["MRC.mx"];
00762 //      }
00763 //      else {
00764                 mrch.mx = nx;
00765 //      }
00766 //      if(dict.has_key("MRC.my")) {
00767 //              mrch.my = dict["MRC.my"];
00768 //      }
00769 //      else {
00770                 mrch.my = ny;
00771 //      }
00772 //      if(dict.has_key("MRC.mz")) {
00773 //              mrch.mz = dict["MRC.mz"];
00774 //      }
00775 //      else {
00776                 mrch.mz = nz;
00777 //      }
00778 
00779         mrch.xlen = mrch.mx * (float) dict["apix_x"];
00780         mrch.ylen = mrch.my * (float) dict["apix_y"];
00781         mrch.zlen = mrch.mz * (float) dict["apix_z"];
00782 
00783         if (dict.has_key("MRC.nxstart")) {
00784                 mrch.nxstart = dict["MRC.nxstart"];
00785         }
00786         else {
00787                 mrch.nxstart = -nx / 2;
00788         }
00789 
00790         if (dict.has_key("MRC.nystart")) {
00791                 mrch.nystart = dict["MRC.nystart"];
00792         }
00793         else {
00794                 mrch.nystart = -ny / 2;
00795         }
00796 
00797         if (dict.has_key("MRC.nzstart")) {
00798                 mrch.nzstart = dict["MRC.nzstart"];
00799         }
00800         else {
00801                 mrch.nzstart = -nz / 2;
00802         }
00803 
00804         strncpy(mrch.map, "MAP ", 4);
00805         mrch.machinestamp = generate_machine_stamp();
00806 
00807         MrcHeader mrch2 = mrch;
00808 
00809         if (opposite_endian || !use_host_endian) {
00810                 swap_header(mrch2);
00811         }
00812 
00813         if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00814                 throw ImageWriteException(filename, "MRC header");
00815         }
00816 
00817         mode_size = get_mode_size(mrch.mode);
00818         is_new_file = false;
00819 
00820         if (is_stack  &&  ! got_one_image) {
00821                 mrch.nz = 1;
00822         }
00823 
00824         // Do not write ctf to mrc header in EMAN2
00825 
00826 //      if( dict.has_key("ctf") ) {
00827 //              vector<float> vctf = dict["ctf"];
00828 //              EMAN1Ctf ctf_;
00829 //              ctf_.from_vector(vctf);
00830 //              write_ctf(ctf_);
00831 //      }
00832 
00833         EXITFUNC;
00834 
00835         return 0;
00836 }
00837 
00838 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool)
00839 {
00840         ENTERFUNC;
00841 
00842         if (! (isFEI || is_stack)) {
00843                 // single image format, index can only be zero
00844 
00845                 image_index = 0;
00846         }
00847 
00848         if (is_transpose && area != 0) {
00849                 printf("Warning: This image dimension is in (y,x,z), "
00850                                 "region I/O not supported, return the whole image instead.");
00851         }
00852 
00853         check_read_access(image_index, rdata);
00854 
00855         if (area && is_complex_mode()) {
00856                 LOGERR("Error: cannot read a region of a complex image.");
00857 
00858                 return 1;
00859         }
00860 
00861         signed char *    scdata = (signed char *)    rdata;
00862         unsigned char *  cdata  = (unsigned char *)  rdata;
00863         short *          sdata  = (short *)          rdata;
00864         unsigned short * usdata = (unsigned short *) rdata;
00865 
00866         size_t size = 0;
00867         int xlen = 0, ylen = 0, zlen = 0;
00868 
00869         if (isFEI) {    // FEI extended MRC
00870                 check_region(area, FloatSize(feimrch.nx, feimrch.ny, feimrch.nz), is_new_file, false);
00871                 portable_fseek(mrcfile, sizeof(MrcHeader)+feimrch.next, SEEK_SET);
00872 
00873                 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00874                                                                   image_index, mode_size,
00875                                                                   feimrch.nx, feimrch.ny, feimrch.nz, area);
00876 
00877                 EMUtil::get_region_dims(area, feimrch.nx, &xlen, feimrch.ny, &ylen, feimrch.nz, &zlen);
00878 
00879                 size = (size_t)xlen * ylen * zlen;
00880         }
00881         else {  // regular MRC
00882                 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00883                 portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00884 
00885                 EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00886                                                                   image_index, mode_size,
00887                                                                   mrch.nx, mrch.ny, mrch.nz, area);
00888 
00889                 EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00890 
00891                 size = (size_t)xlen * ylen * zlen;
00892         }
00893 
00894         if (mrch.mode != MRC_UCHAR  &&  mrch.mode != MRC_CHAR) {
00895                 if (mode_size == sizeof(short)) {
00896                         become_host_endian < short >(sdata, size);
00897                 }
00898                 else if (mode_size == sizeof(float)) {
00899                         become_host_endian < float >(rdata, size);
00900                 }
00901         }
00902 
00903         if (mrch.mode == MRC_UCHAR) {
00904                 for (size_t i = 0; i < size; ++i) {
00905                         size_t j = size - 1 - i;
00906                         // rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
00907                         rdata[j] = static_cast < float >(cdata[j]);
00908                 }
00909         }
00910         else if (mrch.mode == MRC_CHAR) {
00911                 for (size_t i = 0; i < size; ++i) {
00912                         size_t j = size - 1 - i;
00913                         // rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
00914                         rdata[j] = static_cast < float >(scdata[j]);
00915                 }
00916         }
00917         else if (mrch.mode == MRC_SHORT) {
00918                 for (size_t i = 0; i < size; ++i) {
00919                         size_t j = size - 1 - i;
00920                         rdata[j] = static_cast < float >(sdata[j]);
00921                 }
00922         }
00923         else if (mrch.mode == MRC_USHORT) {
00924                 for (size_t i = 0; i < size; ++i) {
00925                         size_t j = size - 1 - i;
00926                         rdata[j] = static_cast < float >(usdata[j]);
00927                 }
00928         }
00929 
00930         if (is_transpose) {
00931                 transpose(rdata, xlen, ylen, zlen);
00932         }
00933 
00934         if (is_complex_mode()) {
00935                 if (! is_ri) {
00936                         Util::ap2ri(rdata, size);
00937                 }
00938 
00939                 Util::flip_complex_phase(rdata, size);
00940                 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00941         }
00942 
00943         EXITFUNC;
00944 
00945         return 0;
00946 }
00947 
00948 int MrcIO::write_data(float *data, int image_index, const Region* area,
00949                                           EMUtil::EMDataType, bool use_host_endian)
00950 {
00951         ENTERFUNC;
00952 
00953         bool append = (image_index == -1);
00954 
00955         if (is_stack  &&  append) {
00956                 image_index = stack_size - 1;
00957         }
00958 
00959         int max_images = 0;
00960 
00961         if (! is_stack) {
00962                 max_images  = 1;
00963                 image_index = 0;
00964         }
00965 
00966         check_write_access(rw_mode, image_index, max_images, data);
00967         check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00968 
00969         int nx, ny, nz;
00970 
00971         if (! area) {
00972                 nx = mrch.nx;
00973                 ny = mrch.ny;
00974                 nz = mrch.nz;
00975         }
00976         else {
00977                 nx = area->get_width();
00978                 ny = area->get_height();
00979                 nz = area->get_depth();
00980         }
00981 
00982         bool got_one_image = (nz > 1);
00983 
00984         if (is_stack  &&  ! got_one_image) {
00985                 nz = 1;
00986         }
00987 
00988         size_t size = (size_t)nx * ny * nz;
00989 
00990         if (is_complex_mode()) {
00991                 nx *= 2;
00992 
00993                 if (! is_ri) {
00994                         Util::ap2ri(data, size);
00995                         is_ri = 1;
00996                 }
00997 
00998                 Util::flip_complex_phase(data, size);
00999                 Util::rotate_phase_origin(data, nx, ny, nz);
01000         }
01001 
01002         portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
01003 
01004         if ((is_big_endian != ByteOrder::is_host_big_endian()) || ! use_host_endian) {
01005                 if (mrch.mode != MRC_UCHAR  &&  mrch.mode != MRC_CHAR) {
01006                         if (mode_size == sizeof(short)) {
01007                                 ByteOrder::swap_bytes((short*) data, size);
01008                         }
01009                         else if (mode_size == sizeof(float)) {
01010                                 ByteOrder::swap_bytes((float*) data, size);
01011                         }
01012                 }
01013         }
01014 
01015         mode_size = get_mode_size(mrch.mode);
01016 
01017 //      int xlen = 0, ylen = 0, zlen = 0;
01018 //      EMUtil::get_region_dims(area, nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
01019 //      int size = xlen * ylen * zlen;
01020 
01021         void * ptr_data = data;
01022 
01023         float rendermin = 0.0f;
01024         float rendermax = 0.0f;
01025 
01026         EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, nz);
01027 
01028         signed char    *  scdata = NULL;
01029         unsigned char  *  cdata  = NULL;
01030         short          *  sdata  = NULL;
01031         unsigned short *  usdata = NULL;
01032 
01033         if (is_stack  &&  ! got_one_image) {
01034                 mrch.nz = stack_size;
01035         }
01036 
01037         if (mrch.mode == MRC_UCHAR) {
01038                 cdata = new unsigned char[size];
01039 
01040                 for (size_t i = 0; i < size; ++i) {
01041                         if (data[i] <= rendermin) {
01042                                 cdata[i] = 0;
01043                         }
01044                         else if (data[i] >= rendermax){
01045                                 cdata[i] = UCHAR_MAX;
01046                         }
01047                         else {
01048                                 cdata[i] = (unsigned char)((data[i] - rendermin) /
01049                                                                 (rendermax - rendermin) *
01050                                                                 UCHAR_MAX);
01051                         }
01052                 }
01053 
01054                 ptr_data = cdata;
01055 
01056                 update_stats(ptr_data, size);
01057         }
01058         else if (mrch.mode == MRC_CHAR) {
01059                 scdata = new signed char[size];
01060 
01061                 for (size_t i = 0; i < size; ++i) {
01062                         if (data[i] <= rendermin) {
01063                                 scdata[i] = SCHAR_MIN;
01064                         }
01065                         else if (data[i] >= rendermax){
01066                                 scdata[i] = SCHAR_MAX;
01067                         }
01068                         else {
01069                                 scdata[i] = (signed char)((data[i] - rendermin) /
01070                                                                 (rendermax - rendermin) *
01071                                                                 (SCHAR_MAX - SCHAR_MIN) + SCHAR_MIN);
01072                         }
01073                 }
01074 
01075                 ptr_data = scdata;
01076 
01077                 update_stats(ptr_data, size);
01078         }
01079         else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
01080                 sdata = new short[size];
01081 
01082                 for (size_t i = 0; i < size; ++i) {
01083                         if (data[i] <= rendermin) {
01084                                 sdata[i] = SHRT_MIN;
01085                         }
01086                         else if (data[i] >= rendermax) {
01087                                 sdata[i] = SHRT_MAX;
01088                         }
01089                         else {
01090                                 sdata[i] = (short)(((data[i] - rendermin) /
01091                                                                 (rendermax - rendermin)) *
01092                                                                 (SHRT_MAX - SHRT_MIN) + SHRT_MIN);
01093                         }
01094                 }
01095 
01096                 ptr_data = sdata;
01097 
01098                 update_stats(ptr_data, size);
01099         }
01100         else if (mrch.mode == MRC_USHORT) {
01101                 usdata = new unsigned short[size];
01102 
01103                 for (size_t i = 0; i < size; ++i) {
01104                         if (data[i] <= rendermin) {
01105                                 usdata[i] = 0;
01106                         }
01107                         else if (data[i] >= rendermax) {
01108                                 usdata[i] = USHRT_MAX;
01109                         }
01110                         else {
01111                                 usdata[i] = (unsigned short)((data[i] - rendermin) /
01112                                                                 (rendermax - rendermin) *
01113                                                                 USHRT_MAX);
01114                         }
01115                 }
01116 
01117                 ptr_data = usdata;
01118 
01119                 update_stats(ptr_data, size);
01120         }
01121 
01122         if (is_stack  &&  ! got_one_image) {
01123                 mrch.nz = 1;
01124         }
01125 
01126         // New way to write data which includes region writing.
01127         // If it is tested to be OK, remove the old code in the
01128         // #if 0  ... #endif block.
01129 
01130         EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
01131                                                           mode_size, mrch.nx, mrch.ny, mrch.nz, area);
01132 
01133         if (cdata)  {delete [] cdata;  cdata  = NULL;}
01134         if (scdata) {delete [] scdata; scdata = NULL;}
01135         if (sdata)  {delete [] sdata;  sdata  = NULL;}
01136         if (usdata) {delete [] usdata; usdata = NULL;}
01137 
01138 #if 0
01139         int row_size = nx * get_mode_size(mrch.mode);
01140         int sec_size = nx * ny;
01141 
01142         unsigned char *  cbuf = new unsigned char[row_size];
01143         unsigned short * sbuf = (unsigned short *) cbuf;
01144 
01145         for (int i = 0; i < nz; i++) {
01146                 int i2 = i * sec_size;
01147 
01148                 for (int j = 0; j < ny; j++) {
01149                         int k = i2 + j * nx;
01150                         void * pbuf = 0;
01151 
01152                         switch (mrch.mode) {
01153                         case MRC_UCHAR:
01154                                 for (int l = 0; l < nx; l++) {
01155                                         cbuf[l] = static_cast < unsigned char >(data[k + l]);
01156                                 }
01157 
01158                                 pbuf = cbuf;
01159                                 fwrite(cbuf, row_size, 1, mrcfile);
01160 
01161                                 break;
01162 
01163                         case MRC_SHORT:
01164                         case MRC_SHORT_COMPLEX:
01165                                 for (int l = 0; l < nx; l++) {
01166                                         sbuf[l] = static_cast < short >(data[k + l]);
01167                                 }
01168 
01169                                 pbuf = sbuf;
01170                                 fwrite(sbuf, row_size, 1, mrcfile);
01171 
01172                                 break;
01173 
01174                         case MRC_USHORT:
01175                                 for (int l = 0; l < nx; l++) {
01176                                         sbuf[l] = static_cast < unsigned short >(data[k + l]);
01177                                 }
01178 
01179                                 pbuf = sbuf;
01180                                 fwrite(sbuf, row_size, 1, mrcfile);
01181 
01182                                 break;
01183 
01184                         case MRC_FLOAT:
01185                         case MRC_FLOAT_COMPLEX:
01186                                 pbuf = &data[k];
01187 
01188                                 break;
01189                         }
01190 
01191                         if (pbuf) {
01192                                 fwrite(pbuf, row_size, 1, mrcfile);
01193                         }
01194                 }
01195         }
01196 
01197         if (cbuf)
01198         {
01199                 delete [] cbuf;
01200                 cbuf = NULL;
01201         }
01202 #endif
01203 
01204         EXITFUNC;
01205         return 0;
01206 }
01207 
01208 void MrcIO::update_stats(void * data, size_t size)
01209 {
01210         float  v;       // variable to hold pixel value
01211         double sum;
01212         double square_sum;
01213         double mean;
01214         double sigma;
01215         double vv;
01216         float  min, max;
01217         
01218         signed char    *  scdata = NULL;
01219         unsigned char  *  cdata  = NULL;
01220         short          *  sdata  = NULL;
01221         unsigned short *  usdata = NULL;
01222 
01223         bool use_schar  = (mrch.mode == MRC_CHAR);
01224         bool use_uchar  = (mrch.mode == MRC_UCHAR);
01225         bool use_short  = (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX);
01226         bool use_ushort = (mrch.mode == MRC_USHORT);
01227         
01228         if (use_uchar) {
01229                 max    = 0.0;
01230                 min    = UCHAR_MAX;
01231                 cdata  = (unsigned char *) data;
01232         }
01233         else if (use_schar) {
01234                 max    = SCHAR_MIN;
01235                 min    = SCHAR_MAX;
01236                 scdata = (signed char *) data;
01237         }
01238         else if (use_short) {
01239                 max    = (float) SHRT_MIN;
01240                 min    = (float) SHRT_MAX;
01241                 sdata  = (short *) data;
01242         }
01243         else if (use_ushort) {
01244                 max    = 0.0f;
01245                 min    = (float) USHRT_MAX;
01246                 usdata = (unsigned short *) data;
01247         }
01248         else {
01249                 throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
01250         }
01251 
01252         sum = 0.0;
01253 
01254         for (size_t i = 0; i < size; i++) {
01255                 if (use_uchar) {
01256                         v = (float) (cdata[i]);
01257                 }
01258                 else if (use_schar) {
01259                         v = (float) (scdata[i]);
01260                 }
01261                 else if (use_short) {
01262                         v = (float) (sdata[i]);
01263                 }
01264                 else {
01265                         v = (float) (usdata[i]);
01266                 }
01267 
01268                 if (v < min) min = v;
01269                 if (v > max) max = v;
01270 
01271                 sum = sum + v;
01272         }
01273 
01274         if (size > 0) {
01275                 mean = sum / (double) size;
01276         }
01277         else {
01278                 mean = 0.0;
01279         }
01280 
01281         square_sum = 0.0;
01282 
01283         for (size_t i = 0; i < size; i++) {
01284                 if (use_uchar) {
01285                         v = (float) (cdata[i]);
01286                 }
01287                 else if (use_schar) {
01288                         v = (float) (scdata[i]);
01289                 }
01290                 else if (use_short) {
01291                         v = (float) (sdata[i]);
01292                 }
01293                 else {
01294                         v = (float) (usdata[i]);
01295                 }
01296 
01297                 vv = v - mean;
01298 
01299                 square_sum = square_sum  +  vv * vv;
01300         }
01301 
01302         if (size > 1) {
01303                 sigma = std::sqrt(square_sum / (double) (size-1));
01304         }
01305         else {
01306                 sigma = 0.0;
01307         }
01308 
01309         /* change mrch.amin / amax / amean / rms here */
01310 
01311         mrch.amin  = min;
01312         mrch.amax  = max;
01313         mrch.amean = mean;
01314         mrch.rms   = sigma;
01315         
01316 //      MrcHeader mrch2 = mrch;
01317 //
01318 // endian issue, can't get use_host_endian argument
01319 //      bool opposite_endian = false;
01320 //
01321 //      if (!is_new_file) {
01322 //              if (is_big_endian != ByteOrder::is_host_big_endian()) {
01323 //                      opposite_endian = true;
01324 //              }
01325 //
01326 //              portable_fseek(mrcfile, 0, SEEK_SET);
01327 //      }
01328 //      
01329 //      if (opposite_endian || !use_host_endian) {
01330 //              swap_header(mrch2);
01331 //      }
01332 
01333         portable_fseek(mrcfile, 0, SEEK_SET);
01334         
01335         if (fwrite(& mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
01336                 throw ImageWriteException(filename, "Error writing MRC header to update statistics.");
01337         }
01338         
01339         portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
01340 }
01341 
01342 bool MrcIO::is_complex_mode()
01343 {
01344         init();
01345 
01346         if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
01347                 return true;
01348         }
01349 
01350         return false;
01351 }
01352 
01353 int MrcIO::read_ctf(Ctf & ctf, int)
01354 {
01355         ENTERFUNC;
01356 
01357         init();
01358 
01359         size_t n = strlen(CTF_MAGIC);
01360         int err = 1;
01361 
01362         if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
01363                 err = ctf.from_string(string(&mrch.labels[0][n]));
01364         }
01365 
01366         EXITFUNC;
01367 
01368         return err;
01369 }
01370 
01371 void MrcIO::write_ctf(const Ctf & ctf, int)
01372 {
01373         ENTERFUNC;
01374 
01375         init();
01376 
01377         string ctf_str = ctf.to_string();
01378 
01379 #ifdef _WIN32
01380         _snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01381 #else
01382          snprintf(&mrch.labels[0][0],80, "%s%s", CTF_MAGIC, ctf_str.c_str());
01383 #endif  //_WIN32
01384 
01385         rewind(mrcfile);
01386 
01387         if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
01388                 throw ImageWriteException(filename, "write CTF info to header failed");
01389         }
01390 
01391         EXITFUNC;
01392 }
01393 
01394 void MrcIO::flush()
01395 {
01396         fflush(mrcfile);
01397 }
01398 
01399 int MrcIO::get_mode_size(int mm)
01400 {
01401         MrcIO::MrcMode m = static_cast < MrcMode > (mm);
01402 
01403         int msize = 0;
01404         switch (m) {
01405         case MRC_CHAR:
01406         case MRC_UCHAR:
01407                 msize = sizeof(char);
01408                 break;
01409         case MRC_SHORT:
01410         case MRC_USHORT:
01411         case MRC_SHORT_COMPLEX:
01412                 msize = sizeof(short);
01413                 break;
01414         case MRC_FLOAT:
01415         case MRC_FLOAT_COMPLEX:
01416                 msize = sizeof(float);
01417                 break;
01418         default:
01419                 msize = 0;
01420         }
01421 
01422         return msize;
01423 }
01424 
01425 int MrcIO::to_em_datatype(int m)
01426 {
01427         EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
01428 
01429         switch (m) {
01430         case MRC_CHAR:
01431                 e = EMUtil::EM_CHAR;
01432                 break;
01433         case MRC_UCHAR:
01434                 e = EMUtil::EM_UCHAR;
01435                 break;
01436         case MRC_SHORT:
01437                 e = EMUtil::EM_SHORT;
01438                 break;
01439         case MRC_USHORT:
01440                 e = EMUtil::EM_USHORT;
01441                 break;
01442         case MRC_SHORT_COMPLEX:
01443                 e = EMUtil::EM_SHORT_COMPLEX;
01444                 break;
01445         case MRC_FLOAT:
01446                 e = EMUtil::EM_FLOAT;
01447                 break;
01448         case MRC_FLOAT_COMPLEX:
01449                 e = EMUtil::EM_FLOAT_COMPLEX;
01450                 break;
01451         default:
01452                 e = EMUtil::EM_UNKNOWN;
01453         }
01454 
01455         return e;
01456 }
01457 
01458 int MrcIO::to_mrcmode(int e, int is_complex)
01459 {
01460         MrcMode m = MRC_UNKNOWN;
01461         EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
01462 
01463         switch (em_type) {
01464         case EMUtil::EM_CHAR:
01465                 m = MRC_CHAR;
01466 
01467                 break;
01468         case EMUtil::EM_UCHAR:
01469                 m = MRC_UCHAR;
01470 
01471                 break;
01472         case EMUtil::EM_USHORT:
01473                 if (is_complex) {
01474                         m = MRC_SHORT_COMPLEX;
01475                 }
01476                 else {
01477                         m = MRC_USHORT;
01478                 }
01479 
01480                 break;
01481         case EMUtil::EM_SHORT:
01482                 if (is_complex) {
01483                         m = MRC_SHORT_COMPLEX;
01484                 }
01485                 else {
01486                         m = MRC_SHORT;
01487                 }
01488 
01489                 break;
01490         case EMUtil::EM_SHORT_COMPLEX:
01491         case EMUtil::EM_USHORT_COMPLEX:
01492                 m = MRC_SHORT_COMPLEX;
01493 
01494                 break;
01495         case EMUtil::EM_INT:
01496         case EMUtil::EM_UINT:
01497         case EMUtil::EM_FLOAT:
01498                 if (is_complex) {
01499                         m = MRC_FLOAT_COMPLEX;
01500                 }
01501                 else {
01502                         m = MRC_FLOAT;
01503                 }
01504 
01505                 break;
01506         case EMUtil::EM_FLOAT_COMPLEX:
01507                 m = MRC_FLOAT_COMPLEX;
01508 
01509                 break;
01510         default:
01511                 m = MRC_FLOAT;
01512         }
01513 
01514         return m;
01515 }
01516 
01517 int MrcIO::generate_machine_stamp()
01518 {
01519         int stamp = 0;
01520         char *p = (char *) (&stamp);
01521 
01522         if (ByteOrder::is_host_big_endian()) {
01523                 p[0] = 0x11;
01524                 p[1] = 0x11;
01525                 p[2] = 0;
01526                 p[3] = 0;
01527         }
01528         else {
01529                 p[0] = 0x44;
01530                 p[1] = 0x41;
01531                 p[2] = 0;
01532                 p[3] = 0;
01533         }
01534 
01535         return stamp;
01536 }
01537 
01538 void MrcIO::swap_header(MrcHeader& mrch)
01539 {
01540         ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
01541         ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
01542 }
01543 
01544 int MrcIO::get_nimg()
01545 {
01546         init();
01547 
01548         return stack_size;
01549 }
01550 
01551 int MrcIO::transpose(float *data, int xlen, int ylen, int zlen) const
01552 {
01553         float * tmp = new float[xlen*ylen];
01554 
01555         for(size_t z=0; z<(size_t)zlen; ++z) {
01556                 for(size_t y=0; y<(size_t)ylen; ++y) {
01557                         for(size_t x=0; x<(size_t)xlen; ++x) {
01558                                 tmp[x*ylen+y] = data[z*xlen*ylen+y*xlen+x];
01559                         }
01560                 }
01561 
01562                 std::copy(tmp, tmp+xlen*ylen, data+z*xlen*ylen);
01563         }
01564 
01565         delete [] tmp;
01566 
01567         return 0;
01568 }