Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members

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 "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                 //become_host_endian((int *) &mrch, NUM_4BYTES_PRE_MAP);
00093                 //become_host_endian((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
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];  //this field specify the extra bytes for symmetry information
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 //#ifndef SPIDERMRC // Spider MRC files don't satisfy the following test
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 //                      return false;
00179                         LOGWARN("image size check fails, still try to read it...");     //when size doesn't match, print error message instead of make it fail
00180                 }
00181                 else {
00182                         return true;
00183                 }
00184 //#endif // SPIDERMRC
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         //single image format, index can only be zero
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; //only read one 2D image from a tilt series
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);        //=1, FEI-MRC file always use short for data type
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;     //offset from end of header to the first dataset
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++) {        //9 tilt angles
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         /* Read extended image header by specified image index*/
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         //remainder 16 4-byte floats not used
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         //single image format, index can only be zero
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 //      if(dict.has_key("MRC.mx")) {
00626 //              mrch.mx = dict["MRC.mx"];
00627 //      }
00628 //      else {
00629                 mrch.mx = nx;
00630 //      }
00631 //      if(dict.has_key("MRC.my")) {
00632 //              mrch.my = dict["MRC.my"];
00633 //      }
00634 //      else {
00635                 mrch.my = ny;
00636 //      }
00637 //      if(dict.has_key("MRC.mz")) {
00638 //              mrch.mz = dict["MRC.mz"];
00639 //      }
00640 //      else {
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                 //single image format, index can only be zero
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) {     //FEI extended MRC, read one 2D image from a tilt series
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 {  //regular MRC
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                         //rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
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         //single image format, index can only be zero
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 //      int xlen = 0, ylen = 0, zlen = 0;
00835 //      EMUtil::get_region_dims(area, nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00836 //      int size = xlen * ylen * zlen;
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         // New way to write data which includes region writing.
00896         // If it is tested to be OK, remove the old code in the
00897         // #if 0  ... #endif block.
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; //variable to hold pixel value
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         /*change mrch.amin/amax/amean.rms here*/
01048         mrch.amin = min;
01049         mrch.amax = max;
01050         mrch.amean = (float)mean;
01051         mrch.rms = sigma;
01052         
01053         MrcHeader mrch2 = mrch;
01054 
01055 //endian issue, can't get use_host_endian argument
01056 //      bool opposite_endian = false;
01057 
01058 //      if (!is_new_file) {
01059 //              if (is_big_endian != ByteOrder::is_host_big_endian()) {
01060 //                      opposite_endian = true;
01061 //              }
01062 //
01063 //              portable_fseek(mrcfile, 0, SEEK_SET);
01064 //      }
01065 //      
01066 //      if (opposite_endian || !use_host_endian) {
01067 //              swap_header(mrch2);
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 }

Generated on Tue Jul 12 13:49:26 2011 for EMAN2 by  doxygen 1.3.9.1