EMAN2
imagicio2.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 #include "imagicio2.h"
00039 #include "portable_fileio.h"
00040 #include "util.h"
00041 #include "geometry.h"
00042 #include "ctf.h"
00043 #include "emassert.h"
00044 #include "transform.h"
00045 
00046 #ifdef _WIN32
00047         #include <ctime>
00048 #endif  //_WIN32
00049 
00050 using namespace EMAN;
00051 
00052 const char *ImagicIO2::HED_EXT = "hed";
00053 const char *ImagicIO2::IMG_EXT = "img";
00054 const char *ImagicIO2::REAL_TYPE_MAGIC = "REAL";
00055 const char *ImagicIO2::CTF_MAGIC = "!-";
00056 
00057 ImagicIO2::ImagicIO2(string file, IOMode rw)
00058 :       filename(file), rw_mode(rw), hed_file(0), img_file(0), initialized(false)
00059 {
00060         hed_filename = Util::change_filename_ext(filename, HED_EXT);
00061         img_filename = Util::change_filename_ext(filename, IMG_EXT);
00062 
00063         is_big_endian = ByteOrder::is_host_big_endian();
00064         is_new_hed = false;
00065         is_new_img = false;
00066         memset(&imagich, 0, sizeof(Imagic4D));
00067         imagich.count = -1;
00068         datatype = IMAGIC_UNKNOWN_TYPE;
00069         nz = 0;
00070 }
00071 
00072 ImagicIO2::~ImagicIO2()
00073 {
00074         if (hed_file) {
00075                 fclose(hed_file);
00076                 hed_file = 0;
00077         }
00078 
00079         if (img_file) {
00080                 fclose(img_file);
00081                 img_file = 0;
00082         }
00083 }
00084 
00085 void ImagicIO2::init()
00086 {
00087         ENTERFUNC;
00088 
00089         if (initialized) {
00090                 return;
00091         }
00092 
00093         initialized = true;
00094 
00095         is_new_hed = false;
00096         is_new_img = false;
00097 
00098         hed_file = sfopen(hed_filename, rw_mode, &is_new_hed);
00099         img_file = sfopen(img_filename, rw_mode, &is_new_img);
00100 
00101         if (is_new_hed != is_new_img) {
00102                 LOGWARN("IMAGIC header file and data file should both exist or both not exist");
00103         }
00104 
00105         if (!is_new_hed) {
00106                 if (fread(&imagich, sizeof(Imagic4D), 1, hed_file) != 1) {
00107                         throw ImageReadException(hed_filename, "IMAGIC4D header");
00108                 }
00109 
00110 //              if (!is_valid(&imagich)) {
00111 //                      throw ImageReadException(hed_filename, "invalid IMAGIC file");
00112 //              }
00113 
00114                 datatype = get_datatype_from_name(imagich.type);
00115 
00116                 if (datatype != IMAGIC_SHORT && datatype != IMAGIC_FLOAT) {
00117                         LOGERR("unsupported imagic data type: %s", imagich.type);
00118                         throw ImageReadException(hed_filename, "unsupported imagic data type");
00119                 }
00120 
00121                 is_big_endian = ByteOrder::is_data_big_endian(&imagich.ny);
00122                 make_header_host_endian(imagich);
00123                 rewind(hed_file);
00124         }
00125 
00126         EXITFUNC;
00127 }
00128 
00129 int ImagicIO2::init_test()
00130 {
00131         ENTERFUNC;
00132 
00133         if (initialized) {
00134                 return 1;
00135         }
00136 
00137         FILE *in = fopen(hed_filename.c_str(), "rb");
00138         if (!in) {
00139                 throw FileAccessException(filename);
00140         }
00141 
00142         char first_block[1024];
00143         size_t n = fread(first_block, sizeof(char), sizeof(first_block), in);
00144 
00145         if (n == 0) {
00146                 LOGERR("file '%s' is an empty file", filename.c_str());
00147                 fclose(in);
00148                 return -1;
00149         }
00150         fclose(in);
00151 
00152         const int *data = reinterpret_cast <const int *>(first_block);
00153         int nx = data[13];
00154         int ny = data[12];
00155         int izold = data[11];
00156 
00157         if(izold==nx*ny) {
00158                 EXITFUNC;
00159                 return -1;      //old style IMAGIC file
00160         }
00161         else {
00162                 EXITFUNC;
00163                 return 0;       //new IMAGIC4D file
00164         }
00165 }
00166 
00167 bool ImagicIO2::is_valid(const void *first_block)
00168 {
00169         ENTERFUNC;
00170 
00171         if (!first_block) {
00172                 return false;
00173         }
00174 
00175         const int *data = static_cast < const int *>(first_block);
00176         int count = data[1];
00177         int headrec = data[3];
00178         int hour = data[7];
00179         int minute = data[8];
00180         int second = data[9];
00181         int rsize = data[10];
00182         int nx = data[13];
00183         int ny = data[12];
00184         int nz = data[60];
00185         int realtype = data[68];
00186 
00187         bool data_big_endian = ByteOrder::is_data_big_endian(&headrec);
00188 
00189         if (data_big_endian != ByteOrder::is_host_big_endian()) {
00190                 ByteOrder::swap_bytes(&count);
00191                 ByteOrder::swap_bytes(&headrec);
00192                 ByteOrder::swap_bytes(&hour);
00193                 ByteOrder::swap_bytes(&rsize);
00194                 ByteOrder::swap_bytes(&nx);
00195                 ByteOrder::swap_bytes(&ny);
00196                 ByteOrder::swap_bytes(&nz);
00197                 ByteOrder::swap_bytes(&realtype);
00198         }
00199 
00200         const int max_dim = 1 << 20;
00201         bool result = false;
00202 
00203         // this field realtype is unique to new Imagic-5 format
00204         if(realtype != VAX_VMS && realtype != LINUX_WINDOWS && realtype != SGI_IBM) {
00205                 EXITFUNC;
00206                 return result;
00207         }
00208         
00209         if (headrec == 1 &&
00210                 count >= 0 && count < max_dim &&
00211                 nx > 0 && nx < max_dim &&
00212                 ny > 0 && ny < max_dim && 
00213                 nz > 0 && nz < max_dim &&
00214                 hour >= 0 && hour < 24 &&
00215                 minute >=0 && minute <60 &&
00216                 second >=0 && second <60) {
00217                 result = true;
00218         }
00219 
00220         EXITFUNC;
00221         return result;
00222 }
00223 
00224 int ImagicIO2::read_header(Dict & dict, int image_index, const Region * area, bool)
00225 {
00226         ENTERFUNC;
00227 
00228         check_read_access(image_index);
00229 
00230         Imagic4D hed;
00231         if (image_index == 0) {
00232                 hed = imagich;
00233         }
00234         else {
00235                 memset(&hed, 0, sizeof(Imagic4D));
00236                 portable_fseek(hed_file, sizeof(Imagic4D) * image_index, SEEK_SET);
00237                 fread(&hed, sizeof(Imagic4D), 1, hed_file);
00238                 make_header_host_endian(hed);
00239         }
00240 
00241         int nz = hed.izlp ? hed.izlp : 1;
00242         check_region(area, FloatSize(hed.nx, hed.ny, nz), is_new_hed, false);
00243 
00244     datatype = get_datatype_from_name(imagich.type);
00245 
00246         int xlen = 0, ylen = 0, zlen = 0;
00247         EMUtil::get_region_dims(area, hed.nx, &xlen, hed.ny, &ylen, nz, &zlen);
00248 
00249         dict["nx"] = xlen;
00250         dict["ny"] = ylen;
00251         dict["nz"] = zlen;
00252 
00253         dict["IMAGIC.imgnum"] = hed.imgnum;
00254         dict["IMAGIC.count"] = hed.count;
00255         dict["IMAGIC.error"] = hed.error;
00256         
00257         dict["IMAGIC.month"] = hed.month;
00258         dict["IMAGIC.day"] = hed.mday;
00259         dict["IMAGIC.year"] = hed.year;
00260         dict["IMAGIC.hour"] = hed.hour;
00261         dict["IMAGIC.minute"] = hed.minute;
00262         dict["IMAGIC.sec"] = hed.sec;
00263 
00264         dict["IMAGIC.rsize"] = hed.rsize;
00265         dict["IMAGIC.izold"] = hed.izold;
00266         
00267         dict["datatype"] = to_em_datatype(datatype);
00268         dict["IMAGIC.type"] = hed.type+'\0';
00269 
00270         dict["IMAGIC.ixold"] = hed.ixold;
00271         dict["IMAGIC.iyold"] = hed.iyold;
00272         
00273         dict["mean"] = hed.avdens;
00274         dict["sigma"] = hed.sigma;
00275 
00276         dict["maximum"] = hed.densmax;
00277         dict["minimum"] = hed.densmin;
00278         
00279         dict["IMAGIC.complex"] = hed.complex;
00280         dict["IMAGIC.defocus1"] = hed.defocus1;
00281         dict["IMAGIC.defocus2"] = hed.defocus2;
00282         dict["IMAGIC.defangle"] = hed.defangle;
00283         dict["IMAGIC.sinostart"] = hed.sinostart;
00284         dict["IMAGIC.sinoend"] = hed.sinoend;
00285 
00286         dict["IMAGIC.label"] = hed.label+'\0';
00287 
00288         dict["IMAGIC.ccc3d"] = hed.ccc3d;
00289         dict["IMAGIC.ref3d"] = hed.ref3d;
00290         dict["IMAGIC.mident"] = hed.mident;
00291         dict["IMAGIC.ezshift"] = hed.ezshift;
00292         
00293         dict["IMAGIC.ealpha"] = hed.ealpha;
00294         dict["IMAGIC.ebeta"] = hed.ebeta;
00295         dict["IMAGIC.egamma"] = hed.egamma;
00296         
00297         dict["IMAGIC.nalisum"] = hed.nalisum;
00298         dict["IMAGIC.pgroup"] = hed.pgroup;
00299         
00300         dict["IMAGIC.i4lp"] = hed.i4lp; //number of objects (1D, 2D, or 3D) in 4D data
00301         
00302         dict["IMAGIC.alpha"] = hed.alpha;
00303         dict["IMAGIC.beta"] = hed.beta;
00304         dict["IMAGIC.gamma"] = hed.gamma;
00305 
00306         dict["IMAGIC.IMAVERS"] = hed.imavers;
00307         dict["IMAGIC.REALTYPE"] = hed.realtype;
00308 
00309         dict["IMAGIC.ANGLE"] = hed.angle;
00310         dict["IMAGIC.VOLTAGE"] = hed.voltage;
00311         dict["IMAGIC.SPABERR"] = hed.spaberr;
00312         dict["IMAGIC.PCOHER"] = hed.pcoher;
00313         dict["IMAGIC.CCC"] = hed.ccc;
00314         dict["IMAGIC.ERRAR"] = hed.errar;
00315         dict["IMAGIC.ERR3D"] = hed.err3d;
00316         dict["IMAGIC.REF"] = hed.ref;
00317         dict["IMAGIC.CLASSNO"] = hed.ref;
00318         dict["IMAGIC.LOCOLD"] = hed.locold;
00319         dict["IMAGIC.REPQUAL"] = hed.repqual;
00320         dict["IMAGIC.ZSHIFT"] = hed.zshift;
00321         dict["IMAGIC.XSHIFT"] = hed.xshift;
00322         dict["IMAGIC.YSHIFT"] = hed.yshift;
00323         dict["IMAGIC.NUMCLS"] = hed.numcls;
00324         dict["IMAGIC.OVQUAL"] = hed.ovqual;
00325         dict["IMAGIC.EANGLE"] = hed.eangle;
00326         dict["IMAGIC.EXSHIFT"] = hed.exshift;
00327         dict["IMAGIC.EYSHIFT"] = hed.eyshift;
00328         dict["IMAGIC.CMTOTVAR"] = hed.cmtotvar;
00329         dict["IMAGIC.INFORMAT"] = hed.informat;
00330         dict["IMAGIC.NUMEIGEN"] = hed.numeigen;
00331         dict["IMAGIC.NIACTIVE"] = hed.niactive;
00332         dict["apix_x"] = hed.resolx;
00333         dict["apix_y"] = hed.resoly;
00334         dict["apix_z"] = hed.resolz;
00335         if(hed.errar==-1.0) {
00336                 dict["IMAGIC.FABOSA1"] = hed.alpha2;
00337                 dict["IMAGIC.FABOSA2"] = hed.beta2;
00338                 dict["IMAGIC.FABOSA3"] = hed.gamma2;
00339         }
00340         else {
00341                 dict["IMAGIC.ALPHA2"] = hed.alpha2;
00342                 dict["IMAGIC.BETA2"] = hed.beta2;
00343                 dict["IMAGIC.GAMMA2"] = hed.gamma2;
00344         }
00345         dict["IMAGIC.NMETRIC"] = hed.nmetric;
00346         dict["IMAGIC.ACTMSA"] = hed.actmsa;
00347 
00348         vector<float> v_coosmsa(hed.coosmsa, hed.coosmsa+69);
00349         dict["IMAGIC.COOSMSA"] = v_coosmsa;
00350 
00351         dict["IMAGIC.EIGVAL"] = hed.coosmsa[19];
00352         dict["IMAGIC.HISTORY"] = hed.history+'\0';
00353 
00354         dict["orientation_convention"] = "IMAGIC";
00355         const float alpha = hed.alpha;
00356         const float beta = hed.beta;
00357         const float gamma = hed.gamma;
00358         dict["euler_alpha"] = alpha;
00359         dict["euler_beta"] = beta;
00360         dict["euler_gamma"] = gamma;
00361         Transform *trans = new Transform();
00362         trans->set_rotation(Dict("type", "imagic", "alpha", alpha, "beta", beta, "gamma", gamma));
00363         if( nz<=1 ) {
00364                 dict["xform.projection"] = trans;
00365         }
00366         else {
00367                 dict["xform.projection"] = trans;
00368                 dict["xform.align3d"] = trans;
00369         }
00370         Ctf * ctf_ = read_ctf(hed);
00371         if( ctf_ != 0) {
00372                 dict["ctf"] = ctf_;
00373         }
00374         if(trans) {delete trans; trans=0;}
00375         if(ctf_) {delete ctf_; ctf_=0;}
00376 
00377         EXITFUNC;
00378         return 0;
00379 }
00380 
00381 ImagicIO2::DataType ImagicIO2::get_datatype_from_name(const char *name) const
00382 {
00383         DataType t = IMAGIC_UNKNOWN_TYPE;
00384 
00385         if (strncmp(name, "PACK",4) == 0) {
00386                 t = IMAGIC_CHAR;
00387         }
00388         else if (strncmp(name, "INTG",4) == 0) {
00389                 t = IMAGIC_SHORT;
00390         }
00391         else if (strncmp(name, REAL_TYPE_MAGIC,4) == 0) {
00392                 t = IMAGIC_FLOAT;
00393         }
00394         else if (strncmp(name, "COMP",4) == 0) {
00395                 t = IMAGIC_FLOAT_COMPLEX;
00396         }
00397         else if (strncmp(name, "RECO",4) == 0) {
00398                 t = IMAGIC_FFT_FLOAT_COMPLEX;
00399         }
00400         return t;
00401 }
00402 
00403 int ImagicIO2::to_em_datatype(DataType t) const
00404 {
00405         switch (t) {
00406         case IMAGIC_CHAR:
00407                 return EMUtil::EM_CHAR;
00408         case IMAGIC_SHORT:
00409                 return EMUtil::EM_SHORT;
00410         case IMAGIC_FLOAT:
00411                 return EMUtil::EM_FLOAT;
00412         case IMAGIC_FLOAT_COMPLEX:
00413         case IMAGIC_FFT_FLOAT_COMPLEX:
00414                 return EMUtil::EM_FLOAT_COMPLEX;
00415         default:
00416                 break;
00417         }
00418 
00419         return EMUtil::EM_UNKNOWN;
00420 }
00421 
00422 Ctf * ImagicIO2::read_ctf(const Imagic4D& hed) const
00423 {
00424         ENTERFUNC;
00425 
00426         Ctf * ctf_ = 0;
00427         size_t n = strlen(CTF_MAGIC);
00428 
00429         if (strncmp(imagich.label, CTF_MAGIC, n) == 0) {
00430                 ctf_ = new EMAN1Ctf();
00431                 string header_label(hed.label);
00432                 // Note: this block was making things crash because it assumed the following if statement
00433                 // was true - I added the if statement (d.woolford)
00434                 if (header_label.size() > 2) {
00435                         string sctf = "O" + header_label.substr(2);
00436                         ctf_->from_string(sctf);
00437                 }
00438         }
00439 
00440         EXITFUNC;
00441         return ctf_;
00442 }
00443 
00444 void ImagicIO2::write_ctf(const Ctf * const ctf, int image_index)
00445 {
00446         ENTERFUNC;
00447         init();
00448 
00449         size_t n = strlen(CTF_MAGIC);
00450         strcpy(imagich.label, CTF_MAGIC);
00451         string ctf_ = ctf->to_string().substr(1);
00452 
00453         //pad ctf_ to 80 char
00454         if(ctf_.size()>80) {
00455                 ctf_ = ctf_.substr(0, 80);
00456         }
00457         else {
00458                 string padded(80 - ctf_.size(), ' ');
00459                 ctf_ = ctf_ + padded;
00460         }
00461 
00462         strncpy(&imagich.label[n], ctf_.c_str(), sizeof(imagich.label) - n);
00463 
00464         rewind(hed_file);
00465         if (fwrite(&imagich, sizeof(Imagic4D), 1, hed_file) != 1) {
00466                 throw ImageWriteException(hed_filename, "Imagic Header");
00467         }
00468 
00469         EXITFUNC;
00470 }
00471 
00472 bool ImagicIO2::is_complex_mode()
00473 {
00474         init();
00475         if (datatype == IMAGIC_FLOAT_COMPLEX || datatype == IMAGIC_FFT_FLOAT_COMPLEX) {
00476                 return true;
00477         }
00478         return false;
00479 }
00480 
00481 void ImagicIO2::flush()
00482 {
00483         fflush(img_file);
00484         fflush(hed_file);
00485 }
00486 
00487 int ImagicIO2::write_header(EMAN::Dict const& dict, int image_index,
00488                 EMAN::Region const* area, EMUtil::EMDataType, bool use_host_endian)
00489 {
00490         ENTERFUNC;
00491 
00492         if(image_index<0) {
00493                 image_index = get_nimg();
00494         }
00495         check_write_access(rw_mode, image_index);
00496 
00497         if (area) {
00498                 check_region(area, FloatSize(imagich.nx, imagich.ny, imagich.izlp),
00499                                          is_new_hed);
00500                 EXITFUNC;
00501                 return 0;
00502         }
00503 
00504         int nx = dict["nx"];
00505         int ny = dict["ny"];
00506         int nz = dict["nz"];
00507         int nimg=0;             //# images currently in file
00508 
00509         if (!is_new_hed) {
00510         datatype = get_datatype_from_name(imagich.type);
00511 
00512                 if (imagich.nx != nx || imagich.ny != ny || imagich.izlp != nz) {
00513                         char desc[256];
00514                         sprintf(desc, "new IMAGIC size %dx%dx%d is not equal to existing size %dx%dx%d",
00515                                         nx, ny, nz, imagich.nx, imagich.ny, imagich.izlp);
00516                         throw ImageWriteException(filename, desc);
00517                 }
00518 
00519         if (datatype!=IMAGIC_FLOAT) {
00520                         throw ImageWriteException(filename, "Attempted write to non REAL Imagic file");
00521                 }
00522 
00523         rewind(hed_file);
00524                 nimg=image_index+1;
00525         }
00526         else {
00527                 nimg = 1;       //new file writing
00528         }
00529 
00530         Imagic4D new_hed;
00531         memset(&new_hed, 0, sizeof(Imagic4D));
00532 
00533         time_t cur_time = time(0);
00534         struct tm *tm = localtime(&cur_time);
00535 
00536         new_hed.error = 0;
00537         new_hed.headrec = 1;    //always 1 a the moment
00538 
00539         new_hed.mday = tm->tm_mday;
00540         new_hed.month = tm->tm_mon;
00541         new_hed.year = tm->tm_year + 1900;
00542         new_hed.hour = tm->tm_hour;
00543         new_hed.minute = tm->tm_min;
00544         new_hed.sec = tm->tm_sec;
00545 
00546         size_t img_size = nx*ny*nz;
00547         if(img_size > (size_t)INT_MAX) {
00548                 new_hed.rsize = -1;
00549         }
00550         else {
00551                 new_hed.rsize = (int)img_size;
00552         }
00553 
00554         new_hed.nx = nx;
00555         new_hed.ny = ny;
00556         new_hed.izlp = nz;
00557 
00558         strncpy(new_hed.type, REAL_TYPE_MAGIC,4);
00559         new_hed.avdens = (float)dict["mean"];
00560         new_hed.sigma = (float)dict["sigma"];
00561         new_hed.densmax = (float)dict["maximum"];
00562         new_hed.densmin = (float)dict["minimum"];
00563 
00564         new_hed.ixold = 0;
00565         new_hed.iyold = 0;
00566 
00567         string new_label = dict.has_key("IMAGIC.label") ? (string) dict["IMAGIC.label"] : string(79, ' ');
00568         sprintf(new_hed.label, "%s", new_label.c_str() );
00569 
00570         string new_history = dict.has_key("IMAGIC.history") ? (string) dict["IMAGIC.history"] : string(227, ' ');
00571         new_history = new_history.substr(0, 228);
00572         sprintf(new_hed.history, "%s", new_history.c_str());
00573 
00574         new_hed.i4lp = nimg;
00575 
00576         Transform * t = 0;
00577         if(nz<=1 && dict.has_key("xform.projection")) {
00578                 t = dict["xform.projection"];
00579         }
00580         else if(nz>1 && dict.has_key("xform.align3d")) {
00581                 t = dict["xform.align3d"];
00582         }
00583 
00584         if(t) {
00585                 Dict d = t->get_rotation("imagic");
00586                 new_hed.alpha = d["alpha"];
00587                 new_hed.beta = d["beta"];
00588                 new_hed.gamma = d["gamma"];
00589                 delete t;
00590                 t=0;
00591         }
00592         else {
00593                 if(dict.has_key("euler_alpha")) new_hed.alpha = dict["euler_alpha"];
00594                 if(dict.has_key("euler_beta")) new_hed.beta = dict["euler_beta"];
00595                 if(dict.has_key("euler_gamma")) new_hed.gamma = dict["euler_gamma"];
00596         }
00597 
00598         new_hed.realtype = generate_machine_stamp();
00599 
00600         new_hed.resolx = dict["apix_x"];
00601         new_hed.resoly = dict["apix_y"];
00602         new_hed.resolz = dict["apix_z"];
00603 
00604         if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian)  swap_header(new_hed);
00605 
00606         // overwrite existing header if necessary
00607         if (image_index>=0 && image_index<nimg) {
00608                 portable_fseek(hed_file, sizeof(Imagic4D)*image_index, SEEK_SET);
00609                 new_hed.imgnum=image_index+1;
00610                 if (is_big_endian != ByteOrder::is_host_big_endian())
00611                                 ByteOrder::swap_bytes((int *) &new_hed.imgnum,1);
00612                 fwrite(&new_hed, sizeof(Imagic4D),1,hed_file);
00613         }
00614 
00615         // update the 1st header with total # images
00616         int ifol = nimg-1;
00617         if (is_big_endian != ByteOrder::is_host_big_endian()) {
00618                 ByteOrder::swap_bytes((int *) &nimg,1);
00619                 ByteOrder::swap_bytes((int *) &ifol,1);
00620         }
00621         portable_fseek(hed_file, sizeof(int), SEEK_SET);
00622         fwrite(&ifol, sizeof(int), 1, hed_file);
00623         portable_fseek(hed_file, 61*sizeof(int), SEEK_SET);     //I4LP(62) is the number of "objects" in file
00624         fwrite(&nimg, sizeof(int), 1, hed_file);
00625 
00626         // header in machine order
00627         if ( (is_big_endian != ByteOrder::is_host_big_endian()) || !use_host_endian)  swap_header(new_hed);
00628         imagich=new_hed;
00629         imagich.count=ifol;
00630         is_new_hed = false;
00631 
00632         if( dict.has_key("ctf") ) {
00633                 Ctf * ctf_ = dict["ctf"];
00634                 write_ctf(ctf_);
00635                 if(ctf_) {delete ctf_; ctf_=0;}
00636         }
00637 
00638         EXITFUNC;
00639         return 0;
00640 }
00641 
00642 int ImagicIO2::generate_machine_stamp() const
00643 {
00644         int machinestamp;
00645 
00646 #ifdef __sgi
00647         machinestamp = SGI_IBM;
00648 #elif defined __OPENVMS__
00649         machinestamp = VAX_VMS;
00650 #else
00651         machinestamp = LINUX_WINDOWS;
00652 #endif
00653 
00654         return machinestamp;
00655 }
00656 
00657 void ImagicIO2::make_header_host_endian(Imagic4D& hed) const
00658 {
00659         if (is_big_endian != ByteOrder::is_host_big_endian()) {
00660                 swap_header(hed);
00661         }
00662 }
00663 
00664 void ImagicIO2::swap_header(Imagic4D & hed) const
00665 {
00666         ByteOrder::swap_bytes((int *) &hed, NUM_4BYTES_PRE_IYLP);
00667         ByteOrder::swap_bytes(&hed.ixold, NUM_4BYTES_AFTER_IXOLD);
00668         ByteOrder::swap_bytes((int *) &hed.ccc3d, NUM_4BYTES_AFTER_NAME);
00669 }
00670 
00671 int ImagicIO2::write_data(float* data, int image_index, const Region * area, EMAN::EMUtil::EMDataType, bool use_host_endian)
00672 {
00673         ENTERFUNC;
00674 
00675         check_write_access(rw_mode, image_index, 0, data);
00676         check_region(area, FloatSize(imagich.nx, imagich.ny, imagich.izlp), is_new_hed);
00677 
00678         if (image_index == -1) {
00679                 portable_fseek(img_file, 0, SEEK_END);
00680         }
00681         else {
00682                 size_t img_size = imagich.nx * imagich.ny * imagich.izlp * sizeof(float);
00683                 portable_fseek(img_file, img_size*image_index, SEEK_SET);
00684         }
00685 
00686         if(is_new_img) {
00687                 if(!use_host_endian) {
00688                         ByteOrder::swap_bytes(data, imagich.nx * imagich.ny * imagich.izlp);
00689                 }
00690         }
00691         else if (is_big_endian != ByteOrder::is_host_big_endian()) {
00692                 ByteOrder::swap_bytes(data, imagich.nx * imagich.ny * imagich.izlp);
00693         }
00694 
00695         EMUtil::process_region_io(data, img_file, WRITE_ONLY, 0,
00696                                                           sizeof(float), imagich.nx, imagich.ny,
00697                                                           imagich.izlp, area, true);
00698 
00699         EXITFUNC;
00700         return 0;
00701 }
00702 
00703 int ImagicIO2::read_data(float* data, int image_index, EMAN::Region const* area, bool)
00704 {
00705         ENTERFUNC;
00706 
00707         check_read_access(image_index, data);
00708         Assert(datatype != IMAGIC_UNKNOWN_TYPE);
00709 
00710         int nx = imagich.ny;
00711         int ny = imagich.nx;
00712         int nz = imagich.izlp ? imagich.izlp : 1;
00713         size_t img_size = (size_t)nx*ny*nz;
00714 
00715         check_region(area, FloatSize(nx, ny, nz), is_new_hed, false);
00716 
00717         portable_fseek(img_file, img_size*image_index*sizeof(float), SEEK_SET);
00718 
00719         short *sdata = (short *) data;
00720         unsigned char *cdata = (unsigned char *) data;
00721         size_t mode_size = get_datatype_size(datatype);
00722 
00725         EMUtil::process_region_io(cdata, img_file, READ_ONLY, 0, mode_size, nx, ny, nz, area, true);
00726 
00727         if (datatype == IMAGIC_FLOAT) {
00728                 become_host_endian(data, img_size);
00729         }
00730         else if (datatype == IMAGIC_SHORT) {
00731                 become_host_endian(sdata, img_size);
00732 
00733                 for (ptrdiff_t j = img_size - 1; j >= 0; --j) {
00734                         data[j] = static_cast < float >(sdata[j]);
00735                 }
00736         }
00737         else {
00738                 throw ImageReadException(filename, "unknown imagic data type");
00739         }
00740 
00741         EXITFUNC;
00742         return 0;
00743 }
00744 
00745 int ImagicIO2::get_nimg()
00746 {
00747         init();
00748         return imagich.count + 1;
00749 }
00750 
00751 bool ImagicIO2::is_image_big_endian()
00752 {
00753         init();
00754         return is_big_endian;
00755 }
00756 
00757 size_t ImagicIO2::get_datatype_size(DataType t) const
00758 {
00759         size_t s = 0;
00760         switch (t) {
00761         case IMAGIC_CHAR:
00762                 s = sizeof(unsigned char);
00763                 break;
00764         case IMAGIC_SHORT:
00765                 s = sizeof(unsigned short);
00766                 break;
00767         case IMAGIC_FLOAT:
00768         case IMAGIC_FLOAT_COMPLEX:
00769         case IMAGIC_FFT_FLOAT_COMPLEX:
00770                 s = sizeof(float);
00771                 break;
00772         default:
00773                 s = 0;
00774         }
00775 
00776         return s;
00777 }