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

hdfio2.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 #ifdef EM_HDF5
00037 
00038 //#define DEBUGHDF      1
00039 
00040 #include "hdfio2.h"
00041 #include "geometry.h"
00042 #include "ctf.h"
00043 #include "emassert.h"
00044 #include "transform.h"
00045 #include "ctf.h"
00046 #include <iostream>
00047 #include <cstring>
00048 
00049 #ifndef WIN32
00050         #include <sys/param.h>
00051 #else
00052         #define  MAXPATHLEN (MAX_PATH * 4)
00053 #endif  //WIN32
00054 
00055 using namespace EMAN;
00056 
00057 static const int ATTR_NAME_LEN = 128;
00058 
00059 HdfIO2::HdfIO2(const string & hdf_filename, IOMode rw)
00060 :       file(-1), group(-1), filename(hdf_filename),
00061         rw_mode(rw), initialized(false)
00062 {
00063         accprop=H5Pcreate(H5P_FILE_ACCESS);
00064 
00065         //STDIO file driver has 2G size limit on 32 bit Linux system
00066         H5Pset_fapl_sec2( accprop );
00067         //H5Pset_fapl_stdio( accprop );
00068 
00069 //      H5Pset_fapl_core( accprop, 1048576, 0  );
00070 //      H5Pset_cache(accprop)
00071         hsize_t dims=1;
00072         simple_space=H5Screate_simple(1,&dims,NULL);
00073 }
00074 
00075 HdfIO2::~HdfIO2()
00076 {
00077         H5Sclose(simple_space);
00078         H5Pclose(accprop);
00079     if (group >= 0) {
00080         H5Gclose(group);
00081     }
00082     if (file >= 0) {
00083                 H5Fflush(file,H5F_SCOPE_GLOBAL);        // If there were no resource leaks, this wouldn't be necessary...
00084                 H5Fclose(file);
00085     }
00086 #ifdef DEBUGHDF
00087         printf("HDf close\n");
00088 #endif
00089 }
00090 
00091 // This reads an already opened attribute and returns the results as an EMObject
00092 // The attribute is not closed
00093 EMObject HdfIO2::read_attr(hid_t attr) {
00094         hid_t type = H5Aget_type(attr);
00095         hid_t spc = H5Aget_space(attr);
00096         H5T_class_t cls = H5Tget_class(type);
00097         size_t sz = H5Tget_size(type);                                          // storage size, arrays handled in the 'space'
00098         hssize_t pts = H5Sget_simple_extent_npoints(spc);       // number of points > 1 if an array of floats or integers
00099 
00100         EMObject ret(0);
00101         int i;
00102 //      unsigned int ui;
00103         float f,*fa;
00104         int * ia;
00105 //      unsigned int * uia;
00106         double d;
00107         char *s;
00108         vector <float> fv((size_t)pts);
00109         vector <int> iv((size_t)pts);
00110 //      vector <unsigned int> uiv(pts);
00111 
00112         float *matrix;
00113         Transform* t;
00114         Ctf* ctf;
00115 //      int r, c, k=0;
00116 
00117         switch (cls) {
00118         case H5T_INTEGER:
00119                 if(pts==1) {
00120                         H5Aread(attr,H5T_NATIVE_INT,&i);
00121                         ret=EMObject(i);
00122                 }
00123                 else {
00124                         ia=(int *)malloc((size_t)pts*sizeof(int));
00125                         H5Aread(attr,H5T_NATIVE_INT,ia);
00126                         for (i=0; i<pts; i++) iv[i]=ia[i];
00127                         free(ia);
00128                         ret=EMObject(iv);
00129                 }
00130                 break;
00131 //      case H5T_UNSIGNED_INTEGER:
00132 //              if(pts==1) {
00133 //                      H5Aread(attr,H5T_NATIVE_UINT,&ui);
00134 //                      ret=EMObject(ui);
00135 //              }
00136 //              else {
00137 //                      uia=(unsigned int *)malloc(pts*sizeof(unsigned int));
00138 //                      H5Aread(attr,H5T_NATIVE_UINT,uia);
00139 //                      for (i=0; i<pts; i++) uiv[i]=uia[i];
00140 //                      free(uia);
00141 //                      ret=EMObject(uiv);
00142 //              }
00143 //              break;
00144         case H5T_FLOAT:
00145                 if (sz==4) {
00146                         if (pts==1) {
00147                                 H5Aread(attr,H5T_NATIVE_FLOAT,&f);
00148                                 ret=EMObject(f);
00149                         }
00150                         else {
00151                                 fa=(float *)malloc((size_t)pts*sizeof(float));
00152                                 H5Aread(attr,H5T_NATIVE_FLOAT,fa);
00153                                 for (i=0; i<pts; i++) fv[i]=fa[i];
00154                                 free(fa);
00155                                 ret=EMObject(fv);
00156                         }
00157                 }
00158                 else if (sz==8) {
00159                         H5Aread(attr,H5T_NATIVE_DOUBLE,&d);
00160                         ret=EMObject(d);
00161                 }
00162                 break;
00163         case H5T_STRING:
00164                 s=(char *)malloc(sz+1);
00165                 H5Aread(attr,type,s);
00166 //              H5Aread(attr,H5T_NATIVE_CHAR,s);
00167                 if(s[0] == 'O' && isdigit(s[1])) {
00168                         ctf = new EMAN1Ctf();
00169                         ctf->from_string(string(s));
00170                         ret = EMObject(ctf);
00171                         delete ctf;
00172                 }
00173                 else if(s[0] == 'E' && isdigit(s[1])) {
00174                         ctf = new EMAN2Ctf();
00175                         ctf->from_string(string(s));
00176                         ret = EMObject(ctf);
00177                         delete ctf;
00178                 }
00179                 else {
00180                         ret=EMObject(s);
00181                 }
00182                 free(s);
00183                 break;
00184         case H5T_COMPOUND:
00185                 matrix = (float*)malloc(12*sizeof(float));
00186                 H5Aread(attr, type, matrix);
00187 //              ret.create_transform3d_by_array(trans3d);
00188                 t = new Transform(matrix);
00189                 ret = EMObject(t);
00190                 free(matrix);
00191                 delete t; t=0;
00192 
00193 //              trans3d = (float*)malloc(16*sizeof(float));     //16 float for a Transform3D object
00194 //              H5Aread(attr, type, trans3d);
00195 // //           ret.create_transform3d_by_array(trans3d);
00196 //              trans = new Transform3D();
00197 //              for(r=0; r<4; ++r) {
00198 //                      for(c=0; c<4; ++c) {
00199 //                              trans->set(r, c, trans3d[k]);
00200 //                              ++k;
00201 //                      }
00202 //              }
00203 //              ret = EMObject(trans);
00204 //              free(trans3d);
00205                 break;
00206         default:
00207                 LOGERR("Unhandled HDF5 metadata %d", cls);
00208         }
00209 
00210         H5Sclose(spc);
00211         H5Tclose(type);
00212 
00213         return ret;
00214 }
00215 
00216 // This writes an attribute with specified name to a given open object
00217 // The attribute is opened and closed. returns 0 on success
00218 int HdfIO2::write_attr(hid_t loc,const char *name,EMObject obj) {
00219         hid_t type=0;
00220         hid_t spc=0;
00221         hsize_t dims=1;
00222         vector <float> fv;
00223         vector <int> iv;
00224         switch(obj.get_type())
00225         {
00226         case EMObject::BOOL: std::cout << "implement this later" << std::endl; break;
00227         case EMObject::INT:
00228                 type=H5Tcopy(H5T_STD_I32LE);
00229                 spc=H5Scopy(simple_space);
00230                 break;
00231         case EMObject::UNSIGNEDINT:
00232                 type=H5Tcopy(H5T_STD_U32LE);
00233                 spc=H5Scopy(simple_space);
00234                 break;
00235         case EMObject::FLOAT:
00236                 type=H5Tcopy(H5T_IEEE_F32LE);
00237                 spc=H5Scopy(simple_space);
00238                 break;
00239         case EMObject::DOUBLE:
00240                 type=H5Tcopy(H5T_IEEE_F64LE);
00241                 spc=H5Scopy(simple_space);
00242                 break;
00243         case EMObject::STRING:
00244         case EMObject::CTF:
00245                 type=H5Tcopy(H5T_C_S1);
00246                 H5Tset_size(type,strlen((const char *)obj)+1);
00247                 spc=H5Screate(H5S_SCALAR);
00248                 break;
00249         case EMObject::FLOATARRAY:
00250                 type=H5Tcopy(H5T_IEEE_F32LE);
00251                 fv=obj;
00252                 dims=fv.size();
00253                 spc=H5Screate_simple(1,&dims,NULL);
00254                 break;
00255         case EMObject::INTARRAY:
00256                 type=H5Tcopy(H5T_STD_I32LE);
00257                 iv=obj;
00258                 dims=iv.size();
00259                 spc=H5Screate_simple(1,&dims,NULL);
00260                 break;
00261 //      case EMObject::TRANSFORM3D:
00262 //              type = H5Tcreate(H5T_COMPOUND, 16 * sizeof(float)); //Transform3D is a 4x4 matrix
00263 //              H5Tinsert(type, "00", 0, H5T_NATIVE_FLOAT);
00264 //              H5Tinsert(type, "01", 1*sizeof(float), H5T_NATIVE_FLOAT);
00265 //              H5Tinsert(type, "02", 2*sizeof(float), H5T_NATIVE_FLOAT);
00266 //              H5Tinsert(type, "03", 3*sizeof(float), H5T_NATIVE_FLOAT);
00267 //              H5Tinsert(type, "10", 4*sizeof(float), H5T_NATIVE_FLOAT);
00268 //              H5Tinsert(type, "11", 5*sizeof(float), H5T_NATIVE_FLOAT);
00269 //              H5Tinsert(type, "12", 6*sizeof(float), H5T_NATIVE_FLOAT);
00270 //              H5Tinsert(type, "13", 7*sizeof(float), H5T_NATIVE_FLOAT);
00271 //              H5Tinsert(type, "20", 8*sizeof(float), H5T_NATIVE_FLOAT);
00272 //              H5Tinsert(type, "21", 9*sizeof(float), H5T_NATIVE_FLOAT);
00273 //              H5Tinsert(type, "22", 10*sizeof(float), H5T_NATIVE_FLOAT);
00274 //              H5Tinsert(type, "23", 11*sizeof(float), H5T_NATIVE_FLOAT);
00275 //              H5Tinsert(type, "30", 12*sizeof(float), H5T_NATIVE_FLOAT);
00276 //              H5Tinsert(type, "31", 13*sizeof(float), H5T_NATIVE_FLOAT);
00277 //              H5Tinsert(type, "32", 14*sizeof(float), H5T_NATIVE_FLOAT);
00278 //              H5Tinsert(type, "33", 15*sizeof(float), H5T_NATIVE_FLOAT);
00279 //              H5Tpack(type);
00280 //
00281 //              dims = 1;       //one compound type
00282 //              spc = H5Screate_simple(1, &dims, NULL);
00283 //              break;
00284         case EMObject::TRANSFORM:
00285                 type = H5Tcreate(H5T_COMPOUND, 12 * sizeof(float)); //Transform is a 3x4 matrix
00286                 H5Tinsert(type, "00", 0, H5T_NATIVE_FLOAT);
00287                 H5Tinsert(type, "01", 1*sizeof(float), H5T_NATIVE_FLOAT);
00288                 H5Tinsert(type, "02", 2*sizeof(float), H5T_NATIVE_FLOAT);
00289                 H5Tinsert(type, "03", 3*sizeof(float), H5T_NATIVE_FLOAT);
00290                 H5Tinsert(type, "10", 4*sizeof(float), H5T_NATIVE_FLOAT);
00291                 H5Tinsert(type, "11", 5*sizeof(float), H5T_NATIVE_FLOAT);
00292                 H5Tinsert(type, "12", 6*sizeof(float), H5T_NATIVE_FLOAT);
00293                 H5Tinsert(type, "13", 7*sizeof(float), H5T_NATIVE_FLOAT);
00294                 H5Tinsert(type, "20", 8*sizeof(float), H5T_NATIVE_FLOAT);
00295                 H5Tinsert(type, "21", 9*sizeof(float), H5T_NATIVE_FLOAT);
00296                 H5Tinsert(type, "22", 10*sizeof(float), H5T_NATIVE_FLOAT);
00297                 H5Tinsert(type, "23", 11*sizeof(float), H5T_NATIVE_FLOAT);
00298                 H5Tpack(type);
00299 
00300                 dims = 1;       //one compound type
00301                 spc = H5Screate_simple(1, &dims, NULL);
00302                 break;
00303         case EMObject::STRINGARRAY:
00304         case EMObject::EMDATA:
00305         case EMObject::XYDATA:
00306         case EMObject::FLOAT_POINTER:
00307         case EMObject::INT_POINTER:
00308         case EMObject::VOID_POINTER:
00309                 return -1;
00310                 break;
00311         case EMObject::UNKNOWN:
00312                 break;
00313         }
00314 
00315     //we need this delete attribute call here, even we called erase_header()
00316     //at the biginning of write_header(), since the  "imageid_max" need be updated correctly.
00317         if( H5Adelete(loc,name) < 0 ) {
00318 #ifdef DEBUGHDF
00319                 LOGERR("Attribute %s deletion error in write_attr().\n", name);
00320 #endif
00321         }
00322         else {
00323 #ifdef DEBUGHDF
00324                 printf("delete attribute %s successfully in write_attr().\n", name);
00325 #endif
00326         }
00327         hid_t attr = H5Acreate(loc,name,type,spc,H5P_DEFAULT);
00328 
00329         unsigned int i;
00330         float f,*fa;
00331         int * ia;
00332         unsigned int ui;
00333         double d;
00334         const char *s;
00335         Transform * tp;
00336         switch(obj.get_type()) {
00337         case EMObject::INT:
00338                 i=(int)obj;
00339                 H5Awrite(attr,H5T_NATIVE_INT,&i);
00340                 break;
00341         case EMObject::UNSIGNEDINT:
00342                 ui=(unsigned int)obj;
00343                 H5Awrite(attr,H5T_NATIVE_UINT,&ui);
00344                 break;
00345         case EMObject::FLOAT:
00346                 f=(float)obj;
00347                 H5Awrite(attr,H5T_NATIVE_FLOAT,&f);
00348                 break;
00349         case EMObject::DOUBLE:
00350                 d=(double)obj;
00351                 H5Awrite(attr,H5T_NATIVE_DOUBLE,&d);
00352                 break;
00353         case EMObject::STRING:
00354         case EMObject::CTF:
00355                 s=(const char *)obj;
00356 //              H5Awrite(attr,H5T_NATIVE_CHAR,s);
00357                 H5Awrite(attr,type,s);
00358                 break;
00359         case EMObject::FLOATARRAY:
00360                 fa=(float *)malloc(fv.size()*sizeof(float));
00361                 for (i=0; i<fv.size(); i++) fa[i]=fv[i];
00362                 H5Awrite(attr,H5T_NATIVE_FLOAT,fa);
00363                 free(fa);
00364                 break;
00365         case EMObject::INTARRAY:
00366                 ia=(int *)malloc(iv.size()*sizeof(int));
00367                 for (i=0; i<iv.size(); i++) ia[i]=iv[i];
00368                 H5Awrite(attr,H5T_NATIVE_INT,ia);
00369                 free(ia);
00370                 break;
00371         case EMObject::TRANSFORM:
00372         {
00373                 tp = (Transform *)obj;
00374                 fa = (float *)malloc(12*sizeof(float));
00375                 int r, c, k=0;
00376                 for(r=0; r<3; ++r) {
00377                         for(c=0; c<4; ++c) {
00378                                 fa[k] = tp->at(r,c);
00379                                 ++k;
00380                         }
00381                 }
00382                 H5Awrite(attr,type,fa);
00383                 free(fa);
00384         }
00385                 break;
00386 //      case EMObject::STRINGARRAY:
00387 //      case EMObject::EMDATA:
00388 //      case EMObject::XYDATA:
00389 //              return -1;
00390 //              break;
00391         default:
00392                 LOGERR("Unhandled HDF5 metadata '%s'", name);
00393 
00394         }
00395 
00396         herr_t ret1 = H5Tclose(type);
00397         herr_t ret2 = H5Sclose(spc);
00398         herr_t ret3 = H5Aclose(attr);
00399         if(ret1>=0 && ret2>=0 && ret3>=0) {
00400                 return 0;
00401         }
00402         else {
00403                 LOGERR("close error in write_attr()\n");
00404                 return -1;
00405         }
00406 }
00407 
00408 // Initializes the file for read-only or read-write access
00409 // Data is stored under /MDF/images
00410 // An attribute named imageid_max stores the number of the highest
00411 // numbered image in the file.
00412 // A group is then made for each individual image, all metadata for the
00413 // individual images is currently associated with the GROUP, not the dataset
00414 // dataset-specific data could also be associated with the dataset in
00415 // future. At the moment, there is only a single dataset in each group.
00416 void HdfIO2::init()
00417 {
00418         ENTERFUNC;
00419         if (initialized) {
00420                 return;
00421         }
00422 #ifdef DEBUGHDF
00423         printf("init\n");
00424 #endif
00425 
00426         H5Eset_auto(0, 0);      // Turn off console error logging.
00427 
00428         if (rw_mode == READ_ONLY) {
00429                 file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, accprop);
00430                 if (file<0) throw FileAccessException(filename);
00431         }
00432         else {
00433                 file = H5Fopen(filename.c_str(), H5F_ACC_RDWR, accprop);
00434                 if (file < 0) {
00435                         file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, accprop);
00436                         if (file < 0) {
00437                                 throw FileAccessException(filename);
00438                         }
00439                         else {
00440 #ifdef DEBUGHDF
00441                                 printf("File truncated or new file created\n");
00442 #endif
00443                         }
00444                 }
00445         }
00446 
00447         group=H5Gopen(file,"/MDF/images");
00448         if (group<0) {
00449                 if (rw_mode == READ_ONLY) throw ImageReadException(filename,"HDF5 file has no image data (no /MDF group)");
00450                 group=H5Gcreate(file,"/MDF",64);                // create the group for Macromolecular data
00451                 if (group<0) throw ImageWriteException(filename,"Unable to add image group (/MDF) to HDF5 file");
00452                 H5Gclose(group);
00453                 group=H5Gcreate(file,"/MDF/images",4096);               // create the group for images/volumes
00454                 if (group<0) throw ImageWriteException(filename,"Unable to add image group (/MDF/images) to HDF5 file");
00455                 write_attr(group,"imageid_max",EMObject(-1));
00456         }
00457         initialized = true;
00458         EXITFUNC;
00459 }
00460 
00461 
00462 // If this version of init() returns -1, then we have an old-style HDF5 file
00463 int HdfIO2::init_test()
00464 {
00465         ENTERFUNC;
00466         if (initialized) {
00467                 return 1;
00468         }
00469 #ifdef DEBUGHDF
00470         printf("init_test\n");
00471 #endif
00472 
00473         H5Eset_auto(0, 0);      // Turn off console error logging.
00474 
00475         hid_t fileid = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5Pcreate(H5P_FILE_ACCESS));
00476         hid_t groupid = H5Gopen(fileid, "/");
00477         hid_t attid = H5Aopen_name(groupid, "num_dataset");
00478 
00479         if (attid < 0) {
00480                 H5Gclose(groupid);
00481                 H5Fclose(fileid);
00482                 init();
00483                 EXITFUNC;
00484                 return 0;
00485         }
00486         else {
00487                 H5Aclose(attid);
00488                 H5Gclose(groupid);
00489                 H5Fclose(fileid);
00490                 EXITFUNC;
00491                 return -1;
00492         }
00493 }
00494 
00495 bool HdfIO2::is_valid(const void *first_block)
00496 {
00497         ENTERFUNC;
00498 
00499         if (first_block) {
00500                 char signature[8] = { 137,72,68,70,13,10,26,10 };
00501                 if (strncmp((const char *)first_block,signature,8)==0) return true;
00502                 // const char* f=(const char *)first_block;
00503                 // printf("bad hdf signature %d %d %d %d %d %d %d %d",f[0],f[1],f[2],f[3],f[4],f[5],f[6],f[7]);
00504                 return false;
00505         }
00506         EXITFUNC;
00507         return false;
00508 }
00509 
00510 // Reads all of the attributes from the /MDF/images/<imgno> group
00511 int HdfIO2::read_header(Dict & dict, int image_index, const Region * area, bool)
00512 {
00513         ENTERFUNC;
00514         init();
00515 #ifdef DEBUGHDF
00516         printf("read_head %d\n", image_index);
00517 #endif
00518         int i;
00519         // Each image is in a group for later expansion. Open the group
00520         char ipath[50];
00521         sprintf(ipath,"/MDF/images/%d", image_index);
00522         hid_t igrp=H5Gopen(file, ipath);
00523 
00524         int nattr=H5Aget_num_attrs(igrp);
00525 
00526         char name[ATTR_NAME_LEN];
00527         for (i=0; i<nattr; i++) {
00528                 hid_t attr=H5Aopen_idx(igrp, i);
00529                 H5Aget_name(attr,127,name);
00530                 if (strncmp(name,"EMAN.", 5)!=0) {
00531                         H5Aclose(attr);
00532                         continue;
00533                 }
00534                 EMObject val=read_attr(attr);
00535                 dict[name+5]=val;
00536                 H5Aclose(attr);
00537         }
00538 
00539         if(dict.has_key("ctf")) {
00540                 string ctfString = (string)dict["ctf"];
00541                 if(ctfString.substr(0, 1) == "O") {
00542                         Ctf * ctf_ = new EMAN1Ctf();
00543                         ctf_->from_string(ctfString);
00544                         dict.erase("ctf");
00545                         dict["ctf"] = ctf_;
00546                         delete ctf_;
00547                 }
00548                 else if(ctfString.substr(0, 1) == "E") {
00549                         Ctf * ctf_ = new EMAN2Ctf();
00550                         ctf_->from_string(ctfString);
00551                         dict.erase("ctf");
00552                         dict["ctf"] = ctf_;
00553                         delete ctf_;
00554                 }
00555         }
00556 
00557         if(area) {
00558                 check_region(area, IntSize(dict["nx"], dict["ny"], dict["nz"]), false, false);
00559 
00560                 dict["nx"] = area->get_width();
00561                 dict["ny"] = area->get_height();
00562                 dict["nz"] = area->get_depth();
00563 
00564                 if( dict.has_key("apix_x") && dict.has_key("apix_y") && dict.has_key("apix_z") )
00565                 {
00566                         if( dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z") )
00567                         {
00568                                 float xorigin = dict["origin_x"];
00569                                 float yorigin = dict["origin_y"];
00570                                 float zorigin = dict["origin_z"];
00571 
00572                                 float apix_x = dict["apix_x"];
00573                                 float apix_y = dict["apix_y"];
00574                                 float apix_z = dict["apix_z"];
00575 
00576                                 dict["origin_x"] = xorigin + apix_x * area->origin[0];
00577                                 dict["origin_y"] = yorigin + apix_y * area->origin[1];
00578                                 dict["origin_z"] = zorigin + apix_z * area->origin[2];
00579                         }
00580                 }
00581         }
00582 
00583         H5Gclose(igrp);
00584         EXITFUNC;
00585         return 0;
00586 }
00587 
00588 // This erases any existing attributes from the image group
00589 // prior to writing a new header. For a new image there
00590 // won't be any, so this should be harmless.
00591 int HdfIO2::erase_header(int image_index)
00592 {
00593         ENTERFUNC;
00594 
00595         if(image_index < 0) return 0; //image_index<0 for appending image, no need for erasing
00596 
00597         init();
00598 #ifdef DEBUGHDF
00599         printf("erase_head %d\n",image_index);
00600 #endif
00601         int i;
00602         // Each image is in a group for later expansion. Open the group
00603         char ipath[50];
00604         sprintf(ipath,"/MDF/images/%d", image_index);
00605         hid_t igrp=H5Gopen(file, ipath);
00606 
00607         int nattr=H5Aget_num_attrs(igrp);
00608 
00609         char name[ATTR_NAME_LEN];
00610         for (i=0; i<nattr; i++) {
00611                 hid_t attr = H5Aopen_idx(igrp, 0); //use 0 as index here, since the H5Adelete() will change the index
00612                 H5Aget_name(attr,127,name);
00613                 H5Aclose(attr);
00614                 if( H5Adelete(igrp,name) < 0 ) {
00615                         LOGERR("attribute %s deletion error in erase_header().\n", name);
00616                 }
00617         }
00618 
00619         H5Gclose(igrp);
00620         EXITFUNC;
00621         return 0;
00622 }
00623 
00624 
00625 int HdfIO2::read_data(float *data, int image_index, const Region *area, bool)
00626 {
00627         ENTERFUNC;
00628 #ifdef DEBUGHDF
00629         printf("read_data %d\n",image_index);
00630 #endif
00631 
00632         char ipath[50];
00633         sprintf(ipath,"/MDF/images/%d/image",image_index);
00634         hid_t ds=H5Dopen(file,ipath);
00635         if (ds<0) throw ImageWriteException(filename,"Image does not exist");
00636         hid_t spc=H5Dget_space(ds);
00637 
00638         hsize_t nx, ny, nz;
00639         hsize_t dims_out[3];
00640         hsize_t rank = H5Sget_simple_extent_ndims(spc);
00641 
00642         H5Sget_simple_extent_dims(spc, dims_out, NULL);
00643         if(rank == 1) {
00644                 nx = dims_out[0];
00645                 ny = 1;
00646                 nz = 1;
00647         }
00648         else if(rank == 2) {
00649                 nx = dims_out[1];
00650                 ny = dims_out[0];
00651                 nz = 1;
00652         }
00653         else if(rank == 3) {
00654                 nx = dims_out[2];
00655                 ny = dims_out[1];
00656                 nz = dims_out[0];
00657         }
00658 
00659         if (area) {
00660                 hid_t memoryspace;
00661 
00662                 /*Get the file dataspace - the region we want to read in the file*/
00663                 int x0, y0, z0;         //the coordinates for up left corner, trim to be within image bound
00664                 int x1, y1, z1;         //the coordinates for down right corner, trim to be within image bound
00665                 int nx1, ny1, nz1;      //dimensions of the sub-region, actual region read form file
00666                 if(rank == 3) {
00667                         hsize_t     doffset[3];             /* hyperslab offset in the file */
00668                         doffset[2] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
00669                         doffset[1] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
00670                         doffset[0] = (hsize_t)(area->z_origin() < 0 ? 0 : area->z_origin());
00671                         x0 = (int)doffset[0];
00672                         y0 = (int)doffset[1];
00673                         z0 = (int)doffset[2];
00674 
00675                         z1 = (int)(area->x_origin() + area->get_width());
00676                         z1 = (int)(z1 > nx ? nx : z1);
00677 
00678                         y1 = (int)(area->y_origin() + area->get_height());
00679                         y1 = (int)(y1 > ny ? ny : y1);
00680                         if(y1 <= 0) {
00681                                 y1 = 1;
00682                         }
00683 
00684                         x1 = (int)(area->z_origin() + area->get_depth());
00685                         x1 = (int)(x1 > nz ? nz : x1);
00686                         if(x1 <= 0) {
00687                                 x1 = 1;
00688                         }
00689 
00690                         if(x1 < x0 || y1< y0 || z1 < z0) return 0; //out of bounds, this is fine, nothing happens
00691 
00692                         hsize_t     dcount[3];              /* size of the hyperslab in the file */
00693                         dcount[0] = x1 - doffset[0];
00694                         dcount[1] = y1 - doffset[1];
00695                         dcount[2] = z1 - doffset[2];
00696 
00697                         H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL, (const hsize_t*)dcount, NULL);
00698 
00699                         /*Define memory dataspace - the memory we will created for the region*/
00700                         hsize_t     dims[3];              /* size of the region in the memory */
00701                         dims[0] = dcount[2]?dcount[2]:1;
00702                         dims[1] = dcount[1]?dcount[1]:1;
00703                         dims[2] = dcount[0]?dcount[0]:1;
00704                         nx1 = (int)dims[0];
00705                         ny1 = (int)dims[1];
00706                         nz1 = (int)dims[2];
00707 
00708                         memoryspace = H5Screate_simple(3, dims, NULL);
00709                 }
00710                 else if(rank == 2) {
00711                         hsize_t     doffset[2];             /* hyperslab offset in the file */
00712                         doffset[1] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
00713                         doffset[0] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
00714                         x0 = (int)doffset[0];
00715                         y0 = (int)doffset[1];
00716                         z0 = 1;
00717 
00718                         y1 = (int)(area->x_origin() + area->get_width());
00719                         y1 = (int)(y1 > nx ? nx : y1);
00720 
00721                         x1 = (int)(area->y_origin() + area->get_height());
00722                         x1 = (int)(x1 > ny ? ny : x1);
00723                         if(x1 <= 0) {
00724                                 x1 = 1;
00725                         }
00726 
00727                         z1 = 1;
00728 
00729                         if(x1 < x0 || y1< y0) return 0; //out of bounds, this is fine, nothing happens
00730 
00731                         hsize_t     dcount[2];              /* size of the hyperslab in the file */
00732                         dcount[0] = x1 - doffset[0];
00733                         dcount[1] = y1 - doffset[1];
00734 
00735                         H5Sselect_none(spc);
00736                         H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL, (const hsize_t*)dcount, NULL);
00737 
00738                         /*Define memory dataspace - the memory we will created for the region*/
00739                         hsize_t     dims[2];              /* size of the region in the memory */
00740                         dims[0] = (hsize_t)(dcount[1]?dcount[1]:1);
00741                         dims[1] = (hsize_t)(dcount[0]?dcount[0]:1);
00742                         nx1 = (int)dims[0];
00743                         ny1 = (int)dims[1];
00744                         nz1 = 1;
00745 
00746                         memoryspace = H5Screate_simple(2, dims, NULL);
00747                 }
00748 
00749                 if( (area->x_origin()>=0) && (area->y_origin()>=0) && (area->z_origin()>=0) && ((hsize_t)(area->x_origin() + area->get_width())<=nx) && ((hsize_t)(area->y_origin() + area->get_height())<=ny) && ((hsize_t)(area->z_origin() + area->get_depth())<=nz) ){      //the region is in boundary
00750                         H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,data);
00751                 }
00752                 else {  //the region are partial out of boundary
00753                         /* When the requested region has some part out of image boundary,
00754                          * we need read the sub-area which is within image,
00755                          * and fill the out of boundary part with zero.
00756                          * We actually read the sub-region from HDF by hyperslab I/O,
00757                          * then copy it back to the pre-allocated region.*/
00758                         float * subdata = new float[nx1*ny1*nz1];
00759 
00760 
00761                         H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,subdata);
00762 
00763                         int xd0=0, yd0=0, zd0=0;        //The coordinates of the top-left corner sub-region in region
00764                         size_t clipped_row_size = 0;
00765                         if(rank == 3) {
00766                                 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
00767                                 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
00768                                 zd0 = (int) (area->z_origin() < 0 ? -area->z_origin() : 0);
00769                                 clipped_row_size = (z1-z0)* sizeof(float);
00770                         }
00771                         else if(rank == 2) {
00772                                 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
00773                                 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
00774                                 clipped_row_size = (y1-y0)* sizeof(float);
00775                         }
00776 
00777                         int src_secsize = nx1 * ny1;
00778                         int dst_secsize = (int)(area->get_width())*(int)(area->get_height());
00779 
00780                         float * src = subdata;
00781                         float * dst = data + zd0*dst_secsize + yd0*(int)(area->get_width()) + xd0;
00782 
00783                         int src_gap = src_secsize - (y1-y0) * nx1;
00784                         int dst_gap = dst_secsize - (y1-y0) * (int)(area->get_width());
00785 
00786                         for(int i = 0; i<nz1; ++i) {
00787                                 for(int j = 0; j<ny1; ++j) {
00788                                         EMUtil::em_memcpy(dst, src, clipped_row_size);
00789 
00790                                         src += nx1;
00791                                         dst += (int)(area->get_width());
00792                                 }
00793                                 src += src_gap;
00794                                 dst += dst_gap;
00795                         }
00796 
00797                         delete [] subdata;
00798                 }
00799                 H5Sclose(memoryspace);
00800         } else {
00801                 H5Dread(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
00802         }
00803 
00804         H5Sclose(spc);
00805         H5Dclose(ds);
00806         EXITFUNC;
00807         return 0;
00808 }
00809 
00810 
00811 // Writes all attributes in 'dict' to the image group
00812 // Creation of the image dataset is also handled here
00813 int HdfIO2::write_header(const Dict & dict, int image_index, const Region* area,
00814                                                 EMUtil::EMDataType, bool)
00815 {
00816 #ifdef DEBUGHDF
00817         printf("write_head %d\n",image_index);
00818 #endif
00819         ENTERFUNC;
00820         init();
00821         
00822         if(image_index<0) {
00823                 image_index = get_nimg();
00824         }
00825 
00826         // If image_index<0 append, and make sure the max value in the file is correct
00827         // though this is normally handled by EMData.write_image()
00828         hid_t attr=H5Aopen_name(group,"imageid_max");
00829         int nimg = read_attr(attr);
00830         H5Aclose(attr);
00831 
00832         unsigned int i;
00833         if (image_index<0) image_index=nimg+1;
00834         if (image_index>nimg) {
00835                 write_attr(group,(const char *)"imageid_max",EMObject(image_index));
00836         }
00837 
00838         // Each image is in a group for later expansion. Open the group, create if necessary
00839         char ipath[50];
00840         sprintf(ipath,"/MDF/images/%d",image_index);
00841         hid_t igrp=H5Gopen(file,ipath);
00842 
00843         if (igrp<0) {   //group not existed
00844                 // Need to create a new image group
00845                 igrp=H5Gcreate(file,ipath,64);          // The image is a group, with attributes on the group
00846                 if (igrp<0) throw ImageWriteException(filename,"Unable to add /MDF/images/# to HDF5 file");
00847 
00848                 sprintf(ipath,"/MDF/images/%d/image",image_index);
00849                 // Now create the actual image dataset
00850                 hid_t space;
00851                 hid_t ds;
00852                 if ((int)dict["nz"]==1)  {
00853                         hsize_t dims[2]= { (int)dict["ny"],(int)dict["nx"] };
00854                         space=H5Screate_simple(2,dims,NULL);
00855                 }
00856                 else {
00857                         hsize_t dims[3]= { (int)dict["nz"],(int)dict["ny"],(int)dict["nx"] };
00858                         space=H5Screate_simple(3,dims,NULL);
00859                 }
00860                 ds=H5Dcreate(file,ipath, H5T_NATIVE_FLOAT, space, H5P_DEFAULT );
00861                 H5Dclose(ds);
00862                 H5Sclose(space);
00863         }
00864         //if group already existed, and the overwiting image is in different size,
00865         //need unlink the existing dataset first
00866         else {
00867                 int nattr=H5Aget_num_attrs(igrp);
00868                 char name[ATTR_NAME_LEN];
00869                 Dict dict2;
00870                 for (int i=0; i<nattr; i++) {
00871                         hid_t attr=H5Aopen_idx(igrp, i);
00872                         H5Aget_name(attr,127,name);
00873                         if (strncmp(name,"EMAN.", 5)!=0) {
00874                                 H5Aclose(attr);
00875                                 continue;
00876                         }
00877                         EMObject val=read_attr(attr);
00878                         dict2[name+5]=val;
00879                         H5Aclose(attr);
00880                 }
00881 
00882                 erase_header(image_index);
00883 
00884                 if((int)dict["nx"]!=(int)dict2["nx"]
00885                                 || (int)dict["ny"]!=(int)dict2["ny"]
00886                                 || (int)dict["nz"]!=(int)dict2["nz"]) {
00887                         sprintf(ipath,"/MDF/images/%d/image",image_index);
00888                         H5Gunlink(igrp, ipath);
00889 
00890                         // Now create the actual image dataset
00891                         hid_t space;
00892                         hid_t ds;
00893                         if ((int)dict["nz"]==1) {
00894                                 hsize_t dims[2]= { (int)dict["ny"],(int)dict["nx"] };
00895                                 space=H5Screate_simple(2,dims,NULL);
00896                         }
00897                         else {
00898                                 hsize_t dims[3]= { (int)dict["nz"],(int)dict["ny"],(int)dict["nx"] };
00899                                 space=H5Screate_simple(3,dims,NULL);
00900                         }
00901                         ds=H5Dcreate(file,ipath, H5T_NATIVE_FLOAT, space, H5P_DEFAULT );
00902                         H5Dclose(ds);
00903                         H5Sclose(space);
00904                 }
00905         }
00906 
00907         if(area) {
00908                 check_region(area, IntSize(dict["nx"], dict["ny"], dict["nz"]), false, true);
00909         }
00910 
00911         // Write the attributes to the group
00912         vector <string> keys=dict.keys();
00913 
00914         for (i=0; i<keys.size(); i++) {
00915                 string s("EMAN.");
00916                 s+=keys[i];
00917                 write_attr(igrp,s.c_str(),dict[keys[i]]);
00918         }
00919 
00920         H5Gclose(igrp);
00921         EXITFUNC;
00922         return 0;
00923 }
00924 
00925 // Writes the actual image data to the corresponding dataset (already created)
00926 int HdfIO2::write_data(float *data, int image_index, const Region* area,
00927                                           EMUtil::EMDataType dt, bool)
00928 {
00929         ENTERFUNC;
00930 
00931 #ifdef DEBUGHDF
00932         printf("write_data %d\n",image_index);
00933 #endif
00934         if (dt!=EMUtil::EM_FLOAT) throw ImageWriteException(filename,"HDF5 write only supports float format");
00935 
00936         if (image_index<0) {
00937                 hid_t attr=H5Aopen_name(group,"imageid_max");
00938                 image_index = read_attr(attr);
00939                 H5Aclose(attr);
00940         }
00941 
00942         char ipath[50];
00943         sprintf(ipath,"/MDF/images/%d/image",image_index);
00944         hid_t ds=H5Dopen(file,ipath);
00945         if (ds<0) throw ImageWriteException(filename,"Image dataset does not exist");
00946 
00947         hid_t spc=H5Dget_space(ds);
00948         if(area) {
00949                 hsize_t doffset[3];             /*hyperslab offset in the file*/
00950                 doffset[0] = (hsize_t)(area->x_origin());
00951                 doffset[1] = (hsize_t)(area->y_origin());
00952                 doffset[2] = (hsize_t)(area->z_origin());
00953 
00954                 hsize_t dcount[3];              /*size of the hyperslab in the file*/
00955                 dcount[0] = (hsize_t)(area->get_width());
00956                 dcount[1] = (hsize_t)(area->get_height()?area->get_height():1);
00957                 dcount[2] = (hsize_t)(area->get_depth()?area->get_depth():1);
00958 
00959                 H5Sselect_hyperslab(spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL, (const hsize_t*)dcount, NULL);
00960 
00961                 /*Create memory space with size of the region.*/
00962                 hsize_t dims[3];        /*size of the region in the memory*/
00963                 dims[0] = (hsize_t)(area->get_width());
00964                 dims[1] = (hsize_t)(area->get_height()?area->get_height():1);
00965                 dims[2] = (hsize_t)(area->get_depth()?area->get_depth():1);
00966 
00967                 hid_t memoryspace = H5Screate_simple(3, dims, NULL);
00968                 H5Dwrite(ds, H5T_NATIVE_FLOAT, memoryspace, spc, H5P_DEFAULT, data);
00969                 H5Sclose(memoryspace);
00970         }
00971         else {
00972                 H5Dwrite(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
00973         }
00974 
00975         H5Sclose(spc);
00976         H5Dclose(ds);
00977         EXITFUNC;
00978         return 0;
00979 }
00980 
00981 int HdfIO2::get_nimg()
00982 {
00983         init();
00984         hid_t attr=H5Aopen_name(group,"imageid_max");
00985         int n = read_attr(attr);
00986         H5Aclose(attr);
00987 
00988         return n+1;
00989 }
00990 
00991 void HdfIO2::flush()
00992 {
00993         return;
00994 }
00995 
00996 bool HdfIO2::is_complex_mode()
00997 {
00998         return false;
00999 }
01000 
01001 // always big endian
01002 bool HdfIO2::is_image_big_endian()
01003 {
01004         return true;
01005 }
01006 
01007 
01008 
01009 #endif  //EM_HDF5

Generated on Fri Apr 30 15:38:51 2010 for EMAN2 by  doxygen 1.3.9.1