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

Generated on Tue Jun 11 12:40:23 2013 for EMAN2 by  doxygen 1.4.7