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

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