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 {
00054         memset(&mrch, 0, sizeof(MrcHeader));
00055         is_ri = 0;
00056         is_big_endian = ByteOrder::is_host_big_endian();
00057         is_new_file = false;
00058         initialized = false;
00059 }
00060 
00061 MrcIO::~MrcIO()
00062 {
00063         if (mrcfile) {
00064                 fclose(mrcfile);
00065                 mrcfile = 0;
00066         }
00067 }
00068 
00069 void MrcIO::init()
00070 {
00071         ENTERFUNC;
00072 
00073         if (initialized) {
00074                 return;
00075         }
00076 
00077         initialized = true;
00078         mrcfile = sfopen(filename, rw_mode, &is_new_file);
00079 
00080         if (!is_new_file) {
00081                 if (fread(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00082                         throw ImageReadException(filename, "MRC header");
00083                 }
00084 
00085                 if (!is_valid(&mrch)) {
00086                         throw ImageReadException(filename, "invalid MRC");
00087                 }
00088 
00089                 is_big_endian = ByteOrder::is_data_big_endian(&mrch.nz);
00090                 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00091                         swap_header(mrch);
00092                 }
00093                 //become_host_endian((int *) &mrch, NUM_4BYTES_PRE_MAP);
00094                 //become_host_endian((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
00095                 mode_size = get_mode_size(mrch.mode);
00096                 if(is_complex_mode()) {
00097                         is_ri = 1;
00098                 }
00099 
00100                 if (mrch.nxstart != 0 || mrch.nystart != 0 || mrch.nzstart != 0) {
00101                         LOGWARN("nx/ny/nz start not zero");
00102                 }
00103 
00104                 if (is_complex_mode()) {
00105                         mrch.nx *= 2;
00106                 }
00107 
00108                 if (mrch.xlen == 0) {
00109                         mrch.xlen = 1.0;
00110                 }
00111 
00112                 if (mrch.ylen == 0) {
00113                         mrch.ylen = 1.0;
00114                 }
00115 
00116                 if (mrch.zlen == 0) {
00117                         mrch.zlen = 1.0;
00118                 }
00119         }
00120         EXITFUNC;
00121 }
00122 
00123 
00124 bool MrcIO::is_image_big_endian()
00125 {
00126         init();
00127         return is_big_endian;
00128 }
00129 
00130 bool MrcIO::is_valid(const void *first_block, off_t file_size)
00131 {
00132         ENTERFUNC;
00133 
00134         if (!first_block) {
00135                 return false;
00136         }
00137 
00138         const int *data = static_cast < const int *>(first_block);
00139         int nx = data[0];
00140         int ny = data[1];
00141         int nz = data[2];
00142         int mrcmode = data[3];
00143         int nsymbt = data[23];  //this field specify the extra bytes for symmetry information
00144 
00145         bool data_big_endian = ByteOrder::is_data_big_endian(&nz);
00146 
00147         if (data_big_endian != ByteOrder::is_host_big_endian()) {
00148                 ByteOrder::swap_bytes(&nx);
00149                 ByteOrder::swap_bytes(&ny);
00150                 ByteOrder::swap_bytes(&nz);
00151                 ByteOrder::swap_bytes(&mrcmode);
00152         }
00153 
00154         if (mrcmode == MRC_SHORT_COMPLEX || mrcmode == MRC_FLOAT_COMPLEX) {
00155                 nx *= 2;
00156         }
00157 
00158         const int max_dim = 1 << 20;
00159 
00160         if ((mrcmode >= MRC_UCHAR && mrcmode < MRC_UNKNOWN) &&
00161                 (nx > 1 && nx < max_dim) && (ny > 0 && ny < max_dim) && (nz > 0 && nz < max_dim)) {
00162 #ifndef SPIDERMRC // Spider MRC files don't satisfy the following test
00163                 if (file_size > 0) {
00164                         off_t file_size1 = (off_t)nx * (off_t)ny * (off_t)nz * (off_t)get_mode_size(mrcmode) + (off_t)sizeof(MrcHeader) + nsymbt;
00165                         if (file_size == file_size1) {
00166                                 return true;
00167                         }
00168                         return false;
00169                 }
00170                 else {
00171                         return true;
00172                 }
00173 #endif // SPIDERMRC
00174                 return true;
00175         }
00176         EXITFUNC;
00177         return false;
00178 }
00179 
00180 int MrcIO::read_header(Dict & dict, int image_index, const Region * area, bool )
00181 {
00182         ENTERFUNC;
00183 
00184         //single image format, index can only be zero
00185         if(image_index == -1) {
00186                 image_index = 0;
00187         }
00188         if(image_index != 0) {
00189                 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().");
00190         }
00191 
00192         init();
00193 
00194         check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file,false);
00195 
00196         dict["apix_x"] = mrch.xlen / mrch.mx;
00197         dict["apix_y"] = mrch.ylen / mrch.my;
00198         dict["apix_z"] = mrch.zlen / mrch.mz;
00199 
00200         dict["MRC.minimum"] = mrch.amin;
00201         dict["MRC.maximum"] = mrch.amax;
00202         dict["MRC.mean"] = mrch.amean;
00203         dict["mean"] = mrch.amean;
00204         dict["datatype"] = to_em_datatype(mrch.mode);
00205 
00206         if (is_complex_mode()) {
00207                 dict["is_complex"] = 1;
00208                 dict["is_complex_ri"] = 1;
00209         }
00210 
00211         int xlen = 0, ylen = 0, zlen = 0;
00212         EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00213 
00214         dict["nx"] = xlen;
00215         dict["ny"] = ylen;
00216         dict["nz"] = zlen;
00217 
00218         if (area) {
00219                 dict["origin_x"] = mrch.xorigin + mrch.xlen * area->origin[0];
00220                 dict["origin_y"] = mrch.yorigin + mrch.xlen * area->origin[1];
00221 
00222                 if (area->get_ndim() == 3 && mrch.nz > 1) {
00223                         dict["origin_z"] = mrch.zorigin + mrch.xlen * area->origin[2];
00224                 }
00225                 else {
00226                         dict["origin_z"] = mrch.zorigin;
00227                 }
00228         }
00229         else {
00230                 dict["origin_x"] = mrch.xorigin;
00231                 dict["origin_y"] = mrch.yorigin;
00232                 dict["origin_z"] = mrch.zorigin;
00233         }
00234 
00235         dict["MRC.nxstart"] = mrch.nxstart;
00236         dict["MRC.nystart"] = mrch.nystart;
00237         dict["MRC.nzstart"] = mrch.nzstart;
00238 
00239         dict["MRC.mx"] = mrch.mx;
00240         dict["MRC.my"] = mrch.my;
00241         dict["MRC.mz"] = mrch.mz;
00242 
00243         dict["MRC.nx"] = mrch.nx;
00244         dict["MRC.ny"] = mrch.ny;
00245         dict["MRC.nz"] = mrch.nz;
00246 
00247         dict["MRC.xlen"] = mrch.xlen;
00248         dict["MRC.ylen"] = mrch.ylen;
00249         dict["MRC.zlen"] = mrch.zlen;
00250 
00251         dict["MRC.alpha"] = mrch.alpha;
00252         dict["MRC.beta"] = mrch.beta;
00253         dict["MRC.gamma"] = mrch.gamma;
00254 
00255         dict["MRC.mapc"] = mrch.mapc;
00256         dict["MRC.mapr"] = mrch.mapr;
00257         dict["MRC.maps"] = mrch.maps;
00258 
00259         dict["MRC.ispg"] = mrch.ispg;
00260         dict["MRC.nsymbt"] = mrch.nsymbt;
00261         dict["MRC.machinestamp"] = mrch.machinestamp;
00262 
00263         dict["MRC.rms"] = mrch.rms;
00264         dict["sigma"] = mrch.rms;
00265         dict["MRC.nlabels"] = mrch.nlabels;
00266         for (int i = 0; i < mrch.nlabels; i++) {
00267                 char label[32];
00268                 sprintf(label, "MRC.label%d", i);
00269                 dict[string(label)] = mrch.labels[i];
00270         }
00271 
00272         EMAN1Ctf ctf_;
00273         if(read_ctf(ctf_) == 0) {
00274                 vector<float> vctf = ctf_.to_vector();
00275                 dict["ctf"] = vctf;
00276         }
00277 
00278         Dict dic;
00279         dic.put("type", "imagic");
00280         dic.put("alpha", mrch.alpha);
00281         dic.put("beta", mrch.beta);
00282         dic.put("gamma", mrch.gamma);
00283         dic.put("tx", mrch.xorigin);
00284         dic.put("ty", mrch.yorigin);
00285         dic.put("tz", mrch.zorigin);
00286         Transform * trans = new Transform(dic);
00287         if(zlen<=1) {
00288                 dict["xform.projection"] = trans;
00289         }
00290         else {
00291                 dict["xform.align3d"] = trans;
00292         }
00293 
00294         if(trans) {delete trans; trans=0;}
00295         EXITFUNC;
00296         return 0;
00297 }
00298 
00299 int MrcIO::write_header(const Dict & dict, int image_index, const Region* area,
00300                                                 EMUtil::EMDataType filestoragetype, bool use_host_endian)
00301 {
00302         ENTERFUNC;
00303 
00304         //single image format, index can only be zero
00305         if(image_index == -1) {
00306                 image_index = 0;
00307         }
00308         if(image_index != 0) {
00309                 throw ImageWriteException(filename, "MRC file does not support stack.");
00310         }
00311         check_write_access(rw_mode, image_index, 1);
00312         if (area) {
00313                 check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00314                 EXITFUNC;
00315                 return 0;
00316         }
00317 
00318         int new_mode = to_mrcmode(filestoragetype, (int) dict["is_complex"]);
00319         int nx = dict["nx"];
00320         int ny = dict["ny"];
00321         int nz = dict["nz"];
00322         is_ri =  dict["is_complex_ri"];
00323 
00324         bool opposite_endian = false;
00325 
00326         if (!is_new_file) {
00327                 if (is_big_endian != ByteOrder::is_host_big_endian()) {
00328                         opposite_endian = true;
00329                 }
00330 #if 0
00331                 if (new_mode != mrch.mode) {
00332                         LOGERR("cannot write to different mode file %s", filename.c_str());
00333                         return 1;
00334                 }
00335 #endif
00336                 portable_fseek(mrcfile, 0, SEEK_SET);
00337         }
00338         else {
00339                 mrch.alpha = mrch.beta = mrch.gamma = 90.0f;
00340                 mrch.mapc = 1;
00341                 mrch.mapr = 2;
00342                 mrch.maps = 3;
00343                 mrch.nxstart = mrch.nystart = mrch.nzstart = 0;
00344         }
00345 
00346         if(nz<=1 && dict.has_key("xform.projection")) {
00347                 Transform * t = dict["xform.projection"];
00348                 Dict d = t->get_params("imagic");
00349                 mrch.alpha = d["alpha"];
00350                 mrch.beta = d["beta"];
00351                 mrch.gamma = d["gamma"];
00352                 mrch.xorigin = d["tx"];
00353                 mrch.yorigin = d["ty"];
00354                 mrch.zorigin = d["tz"];
00355                 if(t) {delete t; t=0;}
00356         }
00357         else if(nz>1 && dict.has_key("xform.align3d")) {
00358                 Transform * t = dict["xform.align3d"];
00359                 Dict d = t->get_params("imagic");
00360                 mrch.alpha = d["alpha"];
00361                 mrch.beta = d["beta"];
00362                 mrch.gamma = d["gamma"];
00363                 mrch.xorigin = d["tx"];
00364                 mrch.yorigin = d["ty"];
00365                 mrch.zorigin = d["tz"];
00366                 if(t) {delete t; t=0;}
00367         }
00368 
00369         if(dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z")){
00370                 mrch.xorigin = (float)dict["origin_x"];
00371                 mrch.yorigin = (float)dict["origin_y"];
00372 
00373                 if (is_new_file) {
00374                         mrch.zorigin = (float)dict["origin_z"];
00375                 }
00376                 else {
00377                         mrch.zorigin = (float) dict["origin_z"] - (float) dict["apix_z"] * image_index;
00378                 }
00379         }
00380 
00381         if (dict.has_key("MRC.nlabels")) {
00382                 mrch.nlabels = dict["MRC.nlabels"];
00383         }
00384 
00385         for (int i = 0; i < MRC_NUM_LABELS; i++) {
00386                 char label[32];
00387                 sprintf(label, "MRC.label%d", i);
00388                 if (dict.has_key(label)) {
00389                         sprintf(&mrch.labels[i][0], "%s", (const char *) dict[label]);
00390                         mrch.nlabels = i + 1;
00391                 }
00392         }
00393 
00394         if (mrch.nlabels < (MRC_NUM_LABELS - 1)) {
00395                 sprintf(&mrch.labels[mrch.nlabels][0], "EMAN %s", Util::get_time_label().c_str());
00396                 mrch.nlabels++;
00397         }
00398 
00399         mrch.labels[mrch.nlabels][0] = '\0';
00400         mrch.mode = new_mode;
00401 
00402         if (is_complex_mode()) {
00403                 mrch.nx = nx / 2;
00404         }
00405         else {
00406                 mrch.nx = nx;
00407         }
00408         mrch.ny = ny;
00409 
00410         if (is_new_file) {
00411                 mrch.nz = nz;
00412         }
00413         else if (image_index >= mrch.nz) {
00414                 mrch.nz = image_index + 1;
00415         }
00416 
00417         mrch.ispg = 0;
00418         mrch.nsymbt = 0;
00419         mrch.amin = dict["minimum"];
00420         mrch.amax = dict["maximum"];
00421         mrch.amean = dict["mean"];
00422         mrch.rms = dict["sigma"];
00423 
00426 //      if(dict.has_key("MRC.mx")) {
00427 //              mrch.mx = dict["MRC.mx"];
00428 //      }
00429 //      else {
00430                 mrch.mx = nx;
00431 //      }
00432 //      if(dict.has_key("MRC.my")) {
00433 //              mrch.my = dict["MRC.my"];
00434 //      }
00435 //      else {
00436                 mrch.my = ny;
00437 //      }
00438 //      if(dict.has_key("MRC.mz")) {
00439 //              mrch.mz = dict["MRC.mz"];
00440 //      }
00441 //      else {
00442                 mrch.mz = nz;
00443 //      }
00444 
00445         mrch.xlen = mrch.mx * (float) dict["apix_x"];
00446         mrch.ylen = mrch.my * (float) dict["apix_y"];
00447         mrch.zlen = mrch.mz * (float) dict["apix_z"];
00448 
00449         if(dict.has_key("MRC.nxstart")) {
00450                 mrch.nxstart = dict["MRC.nxstart"];
00451         }
00452         else {
00453                 mrch.nxstart = -nx / 2;
00454         }
00455         if(dict.has_key("MRC.nystart")) {
00456                 mrch.nystart = dict["MRC.nystart"];
00457         }
00458         else {
00459                 mrch.nystart = -ny / 2;
00460         }
00461         if(dict.has_key("MRC.nzstart")) {
00462                 mrch.nzstart = dict["MRC.nzstart"];
00463         }
00464         else {
00465                 mrch.nzstart = -nz / 2;
00466         }
00467 
00468         sprintf(mrch.map, "MAP ");
00469         mrch.machinestamp = generate_machine_stamp();
00470 
00471         MrcHeader mrch2 = mrch;
00472 
00473         if (opposite_endian || !use_host_endian) {
00474                 swap_header(mrch2);
00475         }
00476 
00477         if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00478                 throw ImageWriteException(filename, "MRC header");
00479         }
00480 
00481         mode_size = get_mode_size(mrch.mode);
00482         is_new_file = false;
00483 
00484         if( dict.has_key("ctf") ) {
00485                 vector<float> vctf = dict["ctf"];
00486                 EMAN1Ctf ctf_;
00487                 ctf_.from_vector(vctf);
00488                 write_ctf(ctf_);
00489         }
00490 
00491         EXITFUNC;
00492         return 0;
00493 }
00494 
00495 int MrcIO::read_data(float *rdata, int image_index, const Region * area, bool )
00496 {
00497         ENTERFUNC;
00498 
00499         //single image format, index can only be zero
00500         image_index = 0;
00501         check_read_access(image_index, rdata);
00502 
00503         if (area && is_complex_mode()) {
00504                 LOGERR("Error: cannot read a region of a complex image.");
00505                 return 1;
00506         }
00507 
00508         check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file, false);
00509 
00510         unsigned char *cdata = (unsigned char *) rdata;
00511         short *sdata = (short *) rdata;
00512         unsigned short *usdata = (unsigned short *) rdata;
00513 
00514         portable_fseek(mrcfile, sizeof(MrcHeader)+mrch.nsymbt, SEEK_SET);
00515 
00516         EMUtil::process_region_io(cdata, mrcfile, READ_ONLY,
00517                                                           image_index, mode_size,
00518                                                           mrch.nx, mrch.ny, mrch.nz, area);
00519 
00520         int xlen = 0, ylen = 0, zlen = 0;
00521         EMUtil::get_region_dims(area, mrch.nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00522 
00523         size_t size = xlen * ylen * zlen;
00524 
00525         if (mrch.mode != MRC_UCHAR) {
00526                 if (mode_size == sizeof(short)) {
00527                         become_host_endian < short >(sdata, size);
00528                 }
00529                 else if (mode_size == sizeof(float)) {
00530                         become_host_endian < float >(rdata, size);
00531                 }
00532         }
00533 
00534         if (mrch.mode == MRC_UCHAR) {
00535                 for (size_t i = 0; i < size; ++i) {
00536                         size_t j = size - 1 - i;
00537                         //rdata[i] = static_cast<float>(cdata[i]/100.0f - 1.28f);
00538                         rdata[j] = static_cast < float >(cdata[j]);
00539                 }
00540         }
00541         else if (mrch.mode == MRC_SHORT ) {
00542                 for (size_t i = 0; i < size; ++i) {
00543                         size_t j = size - 1 - i;
00544                         rdata[j] = static_cast < float >(sdata[j]);
00545                 }
00546         }
00547         else if (mrch.mode == MRC_USHORT) {
00548                 for (size_t i = 0; i < size; ++i) {
00549                         size_t j = size - 1 - i;
00550                         rdata[j] = static_cast < float >(usdata[j]);
00551                 }
00552         }
00553 
00554         if (is_complex_mode()) {
00555                 if(!is_ri) Util::ap2ri(rdata, size);
00556                 Util::flip_complex_phase(rdata, size);
00557                 Util::rotate_phase_origin(rdata, xlen, ylen, zlen);
00558         }
00559         EXITFUNC;
00560         return 0;
00561 }
00562 
00563 int MrcIO::write_data(float *data, int image_index, const Region* area,
00564                                           EMUtil::EMDataType, bool use_host_endian)
00565 {
00566         ENTERFUNC;
00567         //single image format, index can only be zero
00568         image_index = 0;
00569         check_write_access(rw_mode, image_index, 1, data);
00570         check_region(area, FloatSize(mrch.nx, mrch.ny, mrch.nz), is_new_file);
00571 
00572         int nx = mrch.nx;
00573         int ny = mrch.ny;
00574         int nz = mrch.nz;
00575         size_t size = nx * ny * nz;
00576 
00577         if (is_complex_mode()) {
00578                 nx *= 2;
00579                 if (!is_ri) {
00580                         Util::ap2ri(data, size);
00581                         is_ri = 1;
00582                 }
00583                 Util::flip_complex_phase(data, size);
00584                 Util::rotate_phase_origin(data, nx, ny, nz);
00585         }
00586 
00587         portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00588 
00589         if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian) {
00590                 if (mrch.mode != MRC_UCHAR) {
00591                         if (mode_size == sizeof(short)) {
00592                                 ByteOrder::swap_bytes((short*) data, size);
00593                         }
00594                         else if (mode_size == sizeof(float)) {
00595                                 ByteOrder::swap_bytes((float*) data, size);
00596                         }
00597                 }
00598         }
00599         mode_size = get_mode_size(mrch.mode);
00600 
00601 //      int xlen = 0, ylen = 0, zlen = 0;
00602 //      EMUtil::get_region_dims(area, nx, &xlen, mrch.ny, &ylen, mrch.nz, &zlen);
00603 //      int size = xlen * ylen * zlen;
00604         void * ptr_data = data;
00605 
00606         float rendermin = 0.0f;
00607         float rendermax = 0.0f;
00608         EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax);
00609 
00610         unsigned char *cdata = 0;
00611         short *sdata = 0;
00612         unsigned short *usdata = 0;
00613         if (mrch.mode == MRC_UCHAR) {
00614                 cdata = new unsigned char[size];
00615                 for (size_t i = 0; i < size; ++i) {
00616                         if(data[i] <= rendermin) {
00617                                 cdata[i] = 0;
00618                         }
00619                         else if(data[i] >= rendermax){
00620                                 cdata[i] = UCHAR_MAX;
00621                         }
00622                         else {
00623                                 cdata[i]=(unsigned char)((data[i]-rendermin)/(rendermax-rendermin)*UCHAR_MAX);
00624                         }
00625                 }
00626                 ptr_data = cdata;
00627                 update_stat((void *)cdata);
00628         }
00629         else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00630                 sdata = new short[size];
00631                 for (size_t i = 0; i < size; ++i) {
00632                         if(data[i] <= rendermin) {
00633                                 sdata[i] = SHRT_MIN;
00634                         }
00635                         else if(data[i] >= rendermax) {
00636                                 sdata[i] = SHRT_MAX;
00637                         }
00638                         else {
00639                                 sdata[i]=(short)(((data[i]-rendermin)/(rendermax-rendermin))*(SHRT_MAX-SHRT_MIN) - SHRT_MAX);
00640                         }
00641                 }
00642                 ptr_data = sdata;
00643                 update_stat((void *)sdata);
00644         }
00645         else if (mrch.mode == MRC_USHORT) {
00646                 usdata = new unsigned short[size];
00647                 for (size_t i = 0; i < size; ++i) {
00648                         if(data[i] <= rendermin) {
00649                                 usdata[i] = 0;
00650                         }
00651                         else if(data[i] >= rendermax) {
00652                                 usdata[i] = USHRT_MAX;
00653                         }
00654                         else {
00655                                 usdata[i]=(unsigned short)((data[i]-rendermin)/(rendermax-rendermin)*USHRT_MAX);
00656                         }
00657                 }
00658                 ptr_data = usdata;
00659                 update_stat((void *)usdata);
00660         }
00661 
00662         // New way to write data which includes region writing.
00663         // If it is tested to be OK, remove the old code in the
00664         // #if 0  ... #endif block.
00665         EMUtil::process_region_io(ptr_data, mrcfile, WRITE_ONLY, image_index,
00666                                                           mode_size, nx, mrch.ny, mrch.nz, area);
00667 
00668         if(cdata) {delete [] cdata; cdata=0;}
00669         if(sdata) {delete [] sdata; sdata=0;}
00670         if(usdata) {delete [] usdata; usdata=0;}
00671 
00672 #if 0
00673         int row_size = nx * get_mode_size(mrch.mode);
00674         int sec_size = nx * ny;
00675 
00676         unsigned char *cbuf = new unsigned char[row_size];
00677         unsigned short *sbuf = (unsigned short *) cbuf;
00678 
00679         for (int i = 0; i < nz; i++) {
00680                 int i2 = i * sec_size;
00681                 for (int j = 0; j < ny; j++) {
00682                         int k = i2 + j * nx;
00683                         void *pbuf = 0;
00684 
00685                         switch (mrch.mode) {
00686                         case MRC_UCHAR:
00687                                 for (int l = 0; l < nx; l++) {
00688                                         cbuf[l] = static_cast < unsigned char >(data[k + l]);
00689                                 }
00690                                 pbuf = cbuf;
00691                                 fwrite(cbuf, row_size, 1, mrcfile);
00692                                 break;
00693 
00694                         case MRC_SHORT:
00695                         case MRC_SHORT_COMPLEX:
00696                                 for (int l = 0; l < nx; l++) {
00697                                         sbuf[l] = static_cast < short >(data[k + l]);
00698                                 }
00699                                 pbuf = sbuf;
00700                                 fwrite(sbuf, row_size, 1, mrcfile);
00701                                 break;
00702 
00703                         case MRC_USHORT:
00704                                 for (int l = 0; l < nx; l++) {
00705                                         sbuf[l] = static_cast < unsigned short >(data[k + l]);
00706                                 }
00707                                 pbuf = sbuf;
00708                                 fwrite(sbuf, row_size, 1, mrcfile);
00709                                 break;
00710 
00711                         case MRC_FLOAT:
00712                         case MRC_FLOAT_COMPLEX:
00713                                 pbuf = &data[k];
00714                                 break;
00715                         }
00716                         if (pbuf) {
00717                                 fwrite(pbuf, row_size, 1, mrcfile);
00718                         }
00719                 }
00720         }
00721 
00722         if(cbuf)
00723         {
00724                 delete[]cbuf;
00725                 cbuf = 0;
00726         }
00727 #endif
00728 
00729         EXITFUNC;
00730         return 0;
00731 }
00732 
00733 void MrcIO::update_stat(void* data)
00734 {
00735         size_t size =  mrch.nx * mrch.ny * mrch.nz;
00736         float v = 0.0f; //variable to hold pixel value
00737         double sum = 0.0;
00738         double square_sum = 0.0;
00739         double mean = 0.0;
00740         float min, max;
00741         
00742         unsigned char * cdata = 0;
00743         short * sdata = 0;
00744         unsigned short * usdata = 0;
00745         
00746         if (mrch.mode == MRC_UCHAR) {
00747                 max = 0.0f;
00748                 min = UCHAR_MAX;
00749                 cdata = (unsigned char *)data;
00750                 
00751                 for (size_t i = 0; i < size; ++i) {
00752                         v = (float)(cdata[i]);
00753 #ifdef _WIN32
00754                         max = _cpp_max(max,v);
00755                         min = _cpp_min(min,v);
00756 #else
00757                         max=std::max<float>(max,v);
00758                         min=std::min<float>(min,v);
00759 #endif  //_WIN32
00760                         
00761                         sum += v;
00762                         square_sum += v * v;
00763                 }
00764         }
00765         else if (mrch.mode == MRC_SHORT || mrch.mode == MRC_SHORT_COMPLEX) {
00766                 max = (float)SHRT_MIN;
00767                 min = (float)SHRT_MAX;
00768                 sdata = (short *)data;
00769                 
00770                 for (size_t i = 0; i < size; ++i) {
00771                         v = (float)(sdata[i]);
00772 #ifdef _WIN32
00773                         max = _cpp_max(max,v);
00774                         min = _cpp_min(min,v);
00775 #else
00776                         max=std::max<float>(max,v);
00777                         min=std::min<float>(min,v);
00778 #endif  //_WIN32
00779                         
00780                         sum += v;
00781                         square_sum += v * v;
00782                 }
00783         }
00784         else if (mrch.mode == MRC_USHORT) {
00785                 max = 0.0f;
00786                 min = (float)USHRT_MAX;
00787                 usdata = (unsigned short*)data;
00788                 
00789                 for (size_t i = 0; i < size; ++i) {
00790                         v = (float)(usdata[i]);
00791 #ifdef _WIN32
00792                         max = _cpp_max(max,v);
00793                         min = _cpp_min(min,v);
00794 #else
00795                         max=std::max<float>(max,v);
00796                         min=std::min<float>(min,v);
00797 #endif  //_WIN32
00798                         
00799                         sum += v;
00800                         square_sum += v * v;
00801                 }
00802         }
00803         else {
00804                 throw InvalidCallException("This function is used to write 8bit/16bit mrc file only.");
00805         }
00806         
00807         mean = sum/size;
00808 #ifdef _WIN32
00809         float sigma = (float)std::sqrt( _cpp_max(0.0,(square_sum - sum*sum / size)/(size-1)));
00810 #else
00811         float sigma = (float)std::sqrt(std::max<float>(0.0,(square_sum - sum*sum / size)/(size-1)));
00812 #endif  //_WIN32
00813 
00814         /*change mrch.amin/amax/amean.rms here*/
00815         mrch.amin = min;
00816         mrch.amax = max;
00817         mrch.amean = (float)mean;
00818         mrch.rms = sigma;
00819         
00820         MrcHeader mrch2 = mrch;
00821 
00822 //endian issue, can't get use_host_endian argument
00823 //      bool opposite_endian = false;
00824 
00825 //      if (!is_new_file) {
00826 //              if (is_big_endian != ByteOrder::is_host_big_endian()) {
00827 //                      opposite_endian = true;
00828 //              }
00829 //
00830 //              portable_fseek(mrcfile, 0, SEEK_SET);
00831 //      }
00832 //      
00833 //      if (opposite_endian || !use_host_endian) {
00834 //              swap_header(mrch2);
00835 //      }
00836 
00837         portable_fseek(mrcfile, 0, SEEK_SET);
00838         
00839         if (fwrite(&mrch2, sizeof(MrcHeader), 1, mrcfile) != 1) {
00840                 throw ImageWriteException(filename, "MRC header");
00841         }
00842         
00843         portable_fseek(mrcfile, sizeof(MrcHeader), SEEK_SET);
00844 }
00845 
00846 bool MrcIO::is_complex_mode()
00847 {
00848         init();
00849         if (mrch.mode == MRC_SHORT_COMPLEX || mrch.mode == MRC_FLOAT_COMPLEX) {
00850                 return true;
00851         }
00852         return false;
00853 }
00854 
00855 
00856 int MrcIO::read_ctf(Ctf & ctf, int)
00857 {
00858         ENTERFUNC;
00859         init();
00860         size_t n = strlen(CTF_MAGIC);
00861 
00862         int err = 1;
00863         if (strncmp(&mrch.labels[0][0], CTF_MAGIC, n) == 0) {
00864                 err = ctf.from_string(string(&mrch.labels[0][n]));
00865         }
00866         EXITFUNC;
00867         return err;
00868 }
00869 
00870 void MrcIO::write_ctf(const Ctf & ctf, int)
00871 {
00872         ENTERFUNC;
00873         init();
00874 
00875         string ctf_str = ctf.to_string();
00876         sprintf(&mrch.labels[0][0], "%s%s", CTF_MAGIC, ctf_str.c_str());
00877         rewind(mrcfile);
00878 
00879         if (fwrite(&mrch, sizeof(MrcHeader), 1, mrcfile) != 1) {
00880                 throw ImageWriteException(filename, "write CTF info to header failed");
00881         }
00882         EXITFUNC;
00883 }
00884 
00885 void MrcIO::flush()
00886 {
00887         fflush(mrcfile);
00888 }
00889 
00890 
00891 int MrcIO::get_mode_size(int mm)
00892 {
00893         MrcIO::MrcMode m = static_cast < MrcMode > (mm);
00894 
00895         int msize = 0;
00896         switch (m) {
00897         case MRC_UCHAR:
00898                 msize = sizeof(char);
00899                 break;
00900         case MRC_SHORT:
00901         case MRC_USHORT:
00902         case MRC_SHORT_COMPLEX:
00903                 msize = sizeof(short);
00904                 break;
00905         case MRC_FLOAT:
00906         case MRC_FLOAT_COMPLEX:
00907                 msize = sizeof(float);
00908                 break;
00909         default:
00910                 msize = 0;
00911         }
00912 
00913         return msize;
00914 }
00915 
00916 int MrcIO::to_em_datatype(int m)
00917 {
00918         EMUtil::EMDataType e = EMUtil::EM_UNKNOWN;
00919 
00920         switch (m) {
00921         case MRC_UCHAR:
00922                 e = EMUtil::EM_UCHAR;
00923                 break;
00924         case MRC_SHORT:
00925                 e = EMUtil::EM_SHORT;
00926                 break;
00927         case MRC_USHORT:
00928                 e = EMUtil::EM_USHORT;
00929                 break;
00930         case MRC_SHORT_COMPLEX:
00931                 e = EMUtil::EM_SHORT_COMPLEX;
00932                 break;
00933         case MRC_FLOAT:
00934                 e = EMUtil::EM_FLOAT;
00935                 break;
00936         case MRC_FLOAT_COMPLEX:
00937                 e = EMUtil::EM_FLOAT_COMPLEX;
00938                 break;
00939         default:
00940                 e = EMUtil::EM_UNKNOWN;
00941         }
00942         return e;
00943 }
00944 
00945 
00946 int MrcIO::to_mrcmode(int e, int is_complex)
00947 {
00948         MrcMode m = MRC_UNKNOWN;
00949         EMUtil::EMDataType em_type = static_cast < EMUtil::EMDataType > (e);
00950 
00951         switch (em_type) {
00952         case EMUtil::EM_UCHAR:
00953                 m = MRC_UCHAR;
00954                 break;
00955         case EMUtil::EM_USHORT:
00956                 if (is_complex) {
00957                         m = MRC_SHORT_COMPLEX;
00958                 }
00959                 else {
00960                         m = MRC_USHORT;
00961                 }
00962                 break;
00963         case EMUtil::EM_SHORT:
00964                 if (is_complex) {
00965                         m = MRC_SHORT_COMPLEX;
00966                 }
00967                 else {
00968                         m = MRC_SHORT;
00969                 }
00970                 break;
00971         case EMUtil::EM_SHORT_COMPLEX:
00972         case EMUtil::EM_USHORT_COMPLEX:
00973                 m = MRC_SHORT_COMPLEX;
00974                 break;
00975         case EMUtil::EM_CHAR:
00976         case EMUtil::EM_INT:
00977         case EMUtil::EM_UINT:
00978         case EMUtil::EM_FLOAT:
00979                 if (is_complex) {
00980                         m = MRC_FLOAT_COMPLEX;
00981                 }
00982                 else {
00983                         m = MRC_FLOAT;
00984                 }
00985                 break;
00986         case EMUtil::EM_FLOAT_COMPLEX:
00987                 m = MRC_FLOAT_COMPLEX;
00988                 break;
00989         default:
00990                 m = MRC_FLOAT;
00991         }
00992 
00993         return m;
00994 }
00995 
00996 
00997 
00998 int MrcIO::generate_machine_stamp()
00999 {
01000         int stamp = 0;
01001         char *p = (char *) (&stamp);
01002 
01003         if (ByteOrder::is_host_big_endian()) {
01004                 p[0] = 0x11;
01005                 p[1] = 0x11;
01006                 p[2] = 0;
01007                 p[3] = 0;
01008         }
01009         else {
01010                 p[0] = 0x44;
01011                 p[1] = 0x41;
01012                 p[2] = 0;
01013                 p[3] = 0;
01014         }
01015         return stamp;
01016 }
01017 
01018 void MrcIO::swap_header(MrcHeader& mrch)
01019 {
01020         ByteOrder::swap_bytes((int *) &mrch, NUM_4BYTES_PRE_MAP);
01021         ByteOrder::swap_bytes((int *) &mrch.machinestamp, NUM_4BYTES_AFTER_MAP);
01022 }

Generated on Tue May 25 17:34:11 2010 for EMAN2 by  doxygen 1.4.4