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

hdfio.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 #include "hdfio.h"
00039 #include "geometry.h"
00040 #include "ctf.h"
00041 #include "emassert.h"
00042 #include "transform.h"
00043 #include <iostream>
00044 #include <cstring>
00045 
00046 #ifndef _WIN32
00047         #include <sys/param.h>
00048 #else
00049         #include <windows.h>
00050         #define  MAXPATHLEN (MAX_PATH * 4)
00051 #endif  //_WIN32
00052 
00053 using namespace EMAN;
00054 
00055 hid_t HdfIO::mapinfo_type = -1;
00056 const char *HdfIO::HDF5_SIGNATURE = "\211HDF\r\n\032\n";
00057 
00058 herr_t attr_info(hid_t dataset, const char *name, void *opdata)
00059 {
00060         hid_t attr = H5Aopen_name(dataset, name);
00061         float value_float = 0.0f;
00062         int value_int = 0;
00063         string value_string = "";
00064         char * tmp_string = new char[1024];
00065         Dict *dict = (Dict *) opdata;
00066 
00067         if (attr >= 0) {
00068                 hid_t atype = H5Aget_type(attr);
00069                 if (H5Tget_class(atype) == H5T_FLOAT) {
00070                         H5Aread(attr, atype, &value_float);
00071                         (*dict)[name] = value_float;
00072                 }
00073                 else if(H5Tget_class(atype) == H5T_INTEGER) {
00074                         H5Aread(attr, atype, &value_int);
00075                         (*dict)[name] = value_int;
00076                 }
00077                 else if(H5Tget_class(atype) == H5T_STRING) {
00078                         H5Aread(attr, atype, tmp_string);
00079                         value_string = tmp_string;
00080                         (*dict)[name] = value_string;
00081                 }
00082                 else if(H5Tget_class(atype) == H5T_ENUM || H5Tget_class(atype) == H5T_ARRAY) {
00083                         //for old-style hdf file created by EMAN1, the euler convention is enum
00084                         //skip those attribute to make EMAN2 read the content of the image
00085                 }
00086                 else {
00087                         LOGERR("can only handle float/int/string parameters in HDF attr_info()");
00088                         exit(1);
00089                 }
00090                 H5Tclose(atype);
00091                 H5Aclose(attr);
00092         }
00093 
00094         return 0;
00095 }
00096 
00097 HdfIO::HdfIO(const string & hdf_filename, IOMode rw)
00098 :       filename(hdf_filename), rw_mode(rw)
00099 {
00100         initialized = false;
00101         is_new_file = false;
00102         file = -1;
00103         group = -1;
00104         cur_dataset = -1;
00105         cur_image_index = -1;
00106         old_func = 0;
00107         old_client_data = 0;
00108 //      printf("HDf open\n");
00109 }
00110 
00111 HdfIO::~HdfIO()
00112 {
00113         close_cur_dataset();
00114     if (group >= 0) {
00115         H5Gclose(group);
00116     }
00117     if (file >= 0) {
00118                 H5Fflush(file,H5F_SCOPE_GLOBAL);        // If there were no resource leaks, this wouldn't be necessary...
00119                 H5Fclose(file);
00120     }
00121 //printf("HDf close\n");
00122 }
00123 
00124 void HdfIO::init()
00125 {
00126         ENTERFUNC;
00127         if (initialized) {
00128                 return;
00129         }
00130 
00131         initialized = true;
00132 
00133         FILE *tmp_file = sfopen(filename, rw_mode, &is_new_file);
00134 
00135         if (!is_new_file) {
00136                 char buf[128];
00137                 if (fread(buf, sizeof(buf), 1, tmp_file) != 1) {
00138                         fclose(tmp_file);
00139                         throw ImageReadException(filename, "read HDF5 first block");
00140                 }
00141                 else {
00142                         if (!is_valid(buf)) {
00143                                 fclose(tmp_file);
00144                                 throw ImageReadException(filename, "invalid HDF5 file");
00145                         }
00146                 }
00147         }
00148 
00149         fclose(tmp_file);
00150         tmp_file = 0;
00151 
00152         if (rw_mode == READ_ONLY) {
00153                 file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
00154         }
00155         else {
00156                 hdf_err_off();
00157                 file = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
00158                 hdf_err_on();
00159                 if (file < 0) {
00160                         file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);
00161                         H5Fclose(file);
00162                         file = H5Fopen(filename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
00163                 }
00164         }
00165 
00166         if (file < 0) {
00167                 throw FileAccessException(filename);
00168         }
00169 
00170         string root_group_str = get_item_name(ROOT_GROUP);
00171         group = H5Gopen(file, root_group_str.c_str());
00172         cur_dataset = -1;
00173 
00174         H5Giterate(file, root_group_str.c_str(), NULL, file_info, &image_indices);
00175         create_enum_types();
00176         EXITFUNC;
00177 
00178 }
00179 
00180 
00181 bool HdfIO::is_valid(const void *first_block)
00182 {
00183         ENTERFUNC;
00184 
00185         if (first_block) {
00186                 char signature[8] = { 137,72,68,70,13,10,26,10 };
00187                 if (strncmp((const char *)first_block,signature,8)==0) return true;
00188                 // const char* f=(const char *)first_block;
00189                 // 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]);
00190                 return false;
00191         }
00192         EXITFUNC;
00193         return false;
00194 }
00195 
00196 int HdfIO::read_header(Dict & dict, int image_index, const Region * area, bool)
00197 {
00198         ENTERFUNC;
00199 
00200         check_read_access(image_index);
00201 
00202         int nx = 0, ny = 0, nz = 0;
00203         if (get_hdf_dims(image_index, &nx, &ny, &nz) != 0) {
00204                 throw ImageReadException(filename, "invalid image dimensions");
00205         }
00206 
00207         check_region(area, IntSize(nx, ny, nz));
00208 
00209         int xlen = 0, ylen = 0, zlen = 0;
00210         EMUtil::get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
00211 
00212         set_dataset(image_index);
00213         int err = 0;
00214         if (cur_dataset < 0) {
00215                 err = 1;
00216         }
00217         else {
00218                 err = H5Aiterate(cur_dataset, 0, attr_info, &dict);
00219                 if (err < 0) {
00220                         err = 1;
00221                 }
00222         }
00223 
00224         if(dict.has_key("ctf")) {
00225                 Ctf * ctf_ = new EMAN1Ctf();
00226                 ctf_->from_string((string)dict["ctf"]);
00227                 dict.erase("ctf");
00228                 dict["ctf"] = ctf_;
00229                 delete ctf_;
00230         }
00231 
00232         dict["nx"] = xlen;
00233         dict["ny"] = ylen;
00234         dict["nz"] = zlen;
00235 
00236     Dict euler_angles;
00237     read_euler_angles(euler_angles, image_index);
00238 
00239     vector<string> euler_names = euler_angles.keys();
00240     for (size_t i = 0; i < euler_names.size(); i++) {
00241         dict[euler_names[i]] = euler_angles[euler_names[i]];
00242     }
00243 
00244         dict["datatype"] = EMUtil::EM_FLOAT;
00245         EXITFUNC;
00246         return 0;
00247 }
00248 
00249 
00250 int HdfIO::read_data(float *data, int image_index, const Region * area, bool)
00251 {
00252         ENTERFUNC;
00253 
00254         check_read_access(image_index, data);
00255         set_dataset(image_index);
00256 #if 0
00257         if (cur_dataset < 0) {
00258                 char cur_dataset_name[32];
00259                 sprintf(cur_dataset_name, "%d", image_index);
00260                 cur_dataset = H5Dopen(file, cur_dataset_name);
00261 
00262                 if (cur_dataset < 0) {
00263                         char desc[256];
00264                         sprintf(desc, "no image with id = %d", image_index);
00265                         throw ImageReadException(filename, desc);
00266                 }
00267         }
00268 #endif
00269         hid_t datatype = H5Dget_type(cur_dataset);
00270         H5T_class_t t_class = H5Tget_class(datatype);
00271         H5Tclose(datatype);
00272 
00273         if (t_class != H5T_FLOAT) {
00274                 char desc[256];
00275                 sprintf(desc, "unsupported HDF5 data type '%d'", (int) t_class);
00276                 throw ImageReadException(filename, desc);
00277         }
00278 
00279         int nx = 0, ny = 0, nz = 0;
00280         if (get_hdf_dims(image_index, &nx, &ny, &nz) != 0) {
00281                 throw ImageReadException(filename, "invalid image dimensions");
00282         }
00283 
00284         check_region(area, FloatSize(nx, ny, nz));
00285 
00286         int err = 0;
00287         if (!area) {
00288                 err = H5Dread(cur_dataset, H5T_NATIVE_FLOAT, H5S_ALL,
00289                                           H5S_ALL, H5P_DEFAULT, data);
00290         }
00291         else {
00292                 hid_t dataspace_id = 0;
00293                 hid_t memspace_id = 0;
00294 
00295                 err = create_region_space(&dataspace_id, &memspace_id, area,
00296                                                                   nx, ny, nz, image_index);
00297                 if (err == 0) {
00298                         err = H5Dread(cur_dataset, H5T_NATIVE_FLOAT, memspace_id,
00299                                                   dataspace_id, H5P_DEFAULT, data);
00300                 }
00301 
00302                 H5Sclose(dataspace_id);
00303                 H5Sclose(memspace_id);
00304                 if (err < 0) {
00305                         throw ImageReadException(filename,
00306                                                                          "creating memory space or file space id failed");
00307                 }
00308         }
00309 
00310         if (err < 0) {
00311                 char desc[256];
00312                 sprintf(desc, "reading %dth HDF5 image failed", image_index);
00313                 throw ImageReadException(filename, desc);
00314         }
00315 
00316         EXITFUNC;
00317         return 0;
00318 }
00319 
00320 
00321 int HdfIO::write_header(const Dict & dict, int image_index, const Region* area,
00322                                                 EMUtil::EMDataType, bool)
00323 {
00324         ENTERFUNC;
00325         check_write_access(rw_mode, image_index);
00326 
00327         if (area) {
00328                 int nx0 = 0;
00329                 int ny0 = 0;
00330                 int nz0 = 0;
00331 
00332                 if (get_hdf_dims(image_index, &nx0, &ny0, &nz0) != 0) {
00333                         throw ImageReadException(filename, "invalid image dimensions");
00334                 }
00335 
00336                 check_region(area, FloatSize(nx0, ny0, nz0), is_new_file);
00337 
00338                 EXITFUNC;
00339                 return 0;
00340         }
00341 
00342         int nx = dict["nx"];
00343         int ny = dict["ny"];
00344         int nz = dict["nz"];
00345 
00346         create_cur_dataset(image_index, nx, ny, nz);
00347         Assert(cur_dataset >= 0);
00348 
00349         vector<string> keys = dict.keys();
00350         vector<string>::const_iterator iter;
00351         for (iter = keys.begin(); iter != keys.end(); iter++) {
00352                 //handle special case for euler anglers
00353                 if(*iter == "orientation_convention") {
00354                         string eulerstr = (const char*) dict["orientation_convention"];
00355                 write_string_attr(image_index, "orientation_convention", eulerstr);
00356 
00357                 vector<string> euler_names = EMUtil::get_euler_names(eulerstr);
00358                 Dict euler_dict;
00359                 for (size_t i = 0; i < euler_names.size(); i++) {
00360                 euler_dict[euler_names[i]] = dict[euler_names[i]];
00361                 }
00362 
00363                 write_euler_angles(euler_dict, image_index);
00364                 }
00365                 //micrograph_id should save as string
00366                 else if(*iter == "micrograph_id") {
00367                         write_string_attr(image_index, "micrograph_id", (const char *) dict["micrograph_id"]);
00368                 }
00369                 //handle normal attributes
00370                 else {
00371                         EMObject attr_val = dict[*iter];
00372                         //string val_type = EMObject::get_object_type_name(attr_val.get_type());
00373                         EMObject::ObjectType t = attr_val.get_type();
00374                         switch(t) {
00375 
00376                                 case EMObject::INT:
00377                                         write_int_attr(image_index, *iter, attr_val);
00378                                         break;
00379                                 case EMObject::FLOAT:
00380                                 case EMObject::DOUBLE:
00381                                         write_float_attr(image_index, *iter, attr_val);
00382                                         break;
00383                                 case EMObject::STRING:
00384                                         write_string_attr(image_index, *iter, attr_val.to_str());
00385                                         break;
00386 //                              case EMObject::EMDATA:
00387                                 case EMObject::XYDATA:
00388                                 case EMObject::FLOATARRAY:
00389                                 case EMObject::STRINGARRAY:
00390                                         throw NotExistingObjectException("EMObject", "unsupported type");
00391                                         break;
00392                                 case EMObject::UNKNOWN:
00393                                         throw NotExistingObjectException("EMObject", "unsupported type");
00394                                         break;
00395                                 default:
00396                                         throw NotExistingObjectException("EMObject", "unsupported type");
00397                         }
00398                 }
00399         }
00400 
00401 
00402         //flush();
00403         close_cur_dataset();
00404 
00405         EXITFUNC;
00406         return 0;
00407 }
00408 
00409 int HdfIO::write_data(float *data, int image_index, const Region* area,
00410                                           EMUtil::EMDataType, bool)
00411 {
00412         ENTERFUNC;
00413 
00414         check_write_access(rw_mode, image_index, 0, data);
00415 
00416         int nx = read_int_attr(image_index, "nx");
00417         int ny = read_int_attr(image_index, "ny");
00418         int nz = read_int_attr(image_index, "nz");
00419 
00420 
00421         check_region(area, FloatSize(nx, ny, nz), is_new_file);
00422         create_cur_dataset(image_index, nx, ny, nz);
00423         Assert(cur_dataset >= 0);
00424 
00425         int err = 0;
00426 
00427         if (!area) {
00428                 err = H5Dwrite(cur_dataset, H5T_NATIVE_FLOAT, H5S_ALL,
00429                                            H5S_ALL, H5P_DEFAULT, data);
00430                 //if (err >= 0) {
00431                 //increase_num_dataset();
00432                 //image_indices.push_back(image_index);
00433                 //}
00434         }
00435         else {
00436                 hid_t dataspace_id = 0;
00437                 hid_t memspace_id = 0;
00438 
00439                 err = create_region_space(&dataspace_id, &memspace_id, area,
00440                                                                   nx, ny, nz, image_index);
00441                 if (err == 0) {
00442                         err = H5Dwrite(cur_dataset, H5T_NATIVE_FLOAT, memspace_id,
00443                                                    dataspace_id, H5P_DEFAULT, data);
00444                 }
00445 
00446                 H5Sclose(dataspace_id);
00447                 H5Sclose(memspace_id);
00448                 if (err < 0) {
00449                         throw ImageReadException(filename,
00450                                                                          "creating memory space or file space id failed");
00451                 }
00452         }
00453 
00454         if (err < 0) {
00455                 throw ImageWriteException(filename, "HDF data write failed");
00456         }
00457 
00458         //flush();
00459         close_cur_dataset();
00460 
00461         EXITFUNC;
00462         return 0;
00463 }
00464 
00465 void HdfIO::flush()
00466 {
00467         if (cur_dataset > 0) {
00468                 H5Fflush(cur_dataset, H5F_SCOPE_LOCAL);
00469         }
00470 }
00471 
00472 int *HdfIO::read_dims(int image_index, int *p_ndim)
00473 {
00474         set_dataset(image_index);
00475 #if 0
00476         if (cur_dataset < 0) {
00477                 char cur_dataset_name[32];
00478                 sprintf(cur_dataset_name, "%d", image_index);
00479                 cur_dataset = H5Dopen(file, cur_dataset_name);
00480 
00481                 if (cur_dataset < 0) {
00482                         throw ImageReadException(filename, "reading data dimensions");
00483                 }
00484         }
00485 #endif
00486         hid_t dataspace = H5Dget_space(cur_dataset);
00487         int rank = H5Sget_simple_extent_ndims(dataspace);
00488         hsize_t *dims = new hsize_t[rank];
00489         H5Sget_simple_extent_dims(dataspace, dims, NULL);
00490 
00491         int *dims1 = new int[rank];
00492         for (int i = 0; i < rank; i++) {
00493                 dims1[i] = static_cast < int >(dims[i]);
00494         }
00495 
00496         H5Sclose(dataspace);
00497         if( dims )
00498         {
00499                 delete[]dims;
00500                 dims = 0;
00501         }
00502 
00503         (*p_ndim) = rank;
00504         return dims1;
00505 }
00506 
00507 int HdfIO::read_global_int_attr(const string & attr_name)
00508 {
00509         int value = 0;
00510         hid_t attr = H5Aopen_name(group, attr_name.c_str());
00511         if (attr < 0) {
00512                 //LOGWARN("no such hdf attribute '%s'", attr_name.c_str());
00513         }
00514         else {
00515                 H5Aread(attr, H5T_NATIVE_INT, &value);
00516                 H5Aclose(attr);
00517         }
00518 
00519         return value;
00520 }
00521 
00522 float HdfIO::read_global_float_attr(const string & attr_name)
00523 {
00524         float value = 0;
00525         hid_t attr = H5Aopen_name(group, attr_name.c_str());
00526         if (attr >= 0) {
00527                 H5Aread(attr, H5T_NATIVE_INT, &value);
00528                 H5Aclose(attr);
00529         }
00530         return value;
00531 }
00532 
00533 int HdfIO::get_nimg()
00534 {
00535         init();
00536         hdf_err_off();
00537         int n = read_global_int_attr(get_item_name(NUMDATASET));
00538         hdf_err_on();
00539         return n;
00540 }
00541 
00542 
00543 void HdfIO::increase_num_dataset()
00544 {
00545         int n = get_nimg();
00546         n++;
00547         write_global_int_attr(get_item_name(NUMDATASET), n);
00548 }
00549 
00550 
00551 void HdfIO::set_dataset(int image_index)
00552 {
00553         int need_update = 0;
00554 
00555         if (cur_image_index >= 0) {
00556                 if (image_index != cur_image_index) {
00557                         need_update = 1;
00558                 }
00559         }
00560         else {
00561                 cur_image_index = image_index;
00562                 need_update = 1;
00563         }
00564 
00565         if (need_update || cur_dataset < 0) {
00566                 char cur_dataset_name[32];
00567                 sprintf(cur_dataset_name, "%d", image_index);
00568                 hdf_err_off();
00569                 close_cur_dataset();
00570                 cur_dataset = H5Dopen(file, cur_dataset_name);
00571                 if (cur_dataset < 0) {
00572                         throw ImageReadException(filename, "open data set failed");
00573                 }
00574                 hdf_err_on();
00575                 cur_image_index = image_index;
00576         }
00577 
00578         Assert(cur_dataset >= 0);
00579 }
00580 
00581 
00582 
00583 int HdfIO::read_int_attr(int image_index, const string & attr_name)
00584 {
00585         set_dataset(image_index);
00586 
00587         int value = 0;
00588         hid_t attr = H5Aopen_name(cur_dataset, attr_name.c_str());
00589         if (attr >= 0) {
00590                 H5Aread(attr, H5T_NATIVE_INT, &value);
00591                 H5Aclose(attr);
00592         }
00593         else {
00594                 //LOGWARN("no such hdf attribute '%s'", attr_name.c_str());
00595         }
00596 
00597         return value;
00598 }
00599 
00600 
00601 float HdfIO::read_float_attr(int image_index, const string & attr_name)
00602 {
00603         set_dataset(image_index);
00604         return read_float_attr(attr_name);
00605 }
00606 
00607 
00608 float HdfIO::read_float_attr(const string & attr_name)
00609 {
00610         float value = 0;
00611         hid_t attr = H5Aopen_name(cur_dataset, attr_name.c_str());
00612         if (attr >= 0) {
00613                 H5Aread(attr, H5T_NATIVE_FLOAT, &value);
00614                 H5Aclose(attr);
00615         }
00616         else {
00617                 //LOGWARN("no such hdf attribute '%s'", attr_name.c_str());
00618         }
00619 
00620         return value;
00621 }
00622 
00623 
00624 
00625 string HdfIO::read_string_attr(int image_index, const string & attr_name)
00626 {
00627         set_dataset(image_index);
00628         hid_t attr = H5Aopen_name(cur_dataset, attr_name.c_str());
00629 
00630         string value = "";
00631         if (attr < 0) {
00632                 //LOGWARN("no such hdf attribute '%s'", attr_name.c_str());
00633         }
00634         else {
00635                 char *tmp_value = new char[MAXPATHLEN];
00636                 hid_t datatype = H5Tcopy(H5T_C_S1);
00637                 H5Tset_size(datatype, MAXPATHLEN);
00638                 H5Aread(attr, datatype, tmp_value);
00639 
00640                 H5Tclose(datatype);
00641                 H5Aclose(attr);
00642 
00643                 value = tmp_value;
00644 
00645                 if( tmp_value )
00646                 {
00647                         delete[]tmp_value;
00648                         tmp_value = 0;
00649                 }
00650         }
00651 
00652         return value;
00653 }
00654 
00655 int HdfIO::read_array_attr(int image_index, const string & attr_name, void *value)
00656 {
00657         set_dataset(image_index);
00658         int err = 0;
00659 
00660         hid_t attr = H5Aopen_name(cur_dataset, attr_name.c_str());
00661         if (attr < 0) {
00662                 LOGERR("no such hdf attribute '%s'", attr_name.c_str());
00663                 err = 1;
00664         }
00665         else {
00666                 hid_t datatype = H5Aget_type(attr);
00667                 H5Aread(attr, datatype, value);
00668                 H5Tclose(datatype);
00669                 H5Aclose(attr);
00670         }
00671         return err;
00672 }
00673 #if 0
00674 float HdfIO::read_euler_attr(int image_index, const string &attr_name)
00675 {
00676     string full_attr_name = string(EULER_MAGIC) + attr_name;
00677     return read_float_attr(image_index, full_attr_name);
00678 }
00679 #endif
00680 
00681 int HdfIO::read_mapinfo_attr(int image_index, const string & attr_name)
00682 {
00683         set_dataset(image_index);
00684 
00685         MapInfoType val = ICOS_UNKNOWN;
00686         hid_t attr = H5Aopen_name(cur_dataset, attr_name.c_str());
00687         if (attr < 0) {
00688                 LOGERR("no such hdf attribute '%s'", attr_name.c_str());
00689         }
00690         else {
00691                 H5Aread(attr, mapinfo_type, &val);
00692                 H5Aclose(attr);
00693         }
00694         return static_cast < int >(val);
00695 }
00696 
00697 int HdfIO::write_int_attr(int image_index, const string & attr_name, int value)
00698 {
00699         set_dataset(image_index);
00700         return write_int_attr(attr_name, value);
00701 }
00702 
00703 
00704 int HdfIO::write_int_attr(const string & attr_name, int value)
00705 {
00706         int err = -1;
00707         delete_attr(attr_name);
00708 
00709         hid_t dataspace = H5Screate(H5S_SCALAR);
00710         hid_t attr = H5Acreate(cur_dataset, attr_name.c_str(), H5T_NATIVE_INT, dataspace, H5P_DEFAULT);
00711 
00712         if (attr >= 0) {
00713                 err = H5Awrite(attr, H5T_NATIVE_INT, &value);
00714         }
00715 
00716         H5Aclose(attr);
00717         H5Sclose(dataspace);
00718 
00719         if (err < 0) {
00720                 return 1;
00721         }
00722 
00723         return 0;
00724 }
00725 
00726 int HdfIO::write_float_attr_from_dict(int image_index, const string & attr_name,
00727                                                                           const Dict & dict)
00728 {
00729         if (dict.has_key(attr_name)) {
00730                 return write_float_attr(image_index, attr_name, dict[attr_name]);
00731         }
00732         return 0;
00733 }
00734 
00735 int HdfIO::write_float_attr(int image_index, const string & attr_name, float value)
00736 {
00737         set_dataset(image_index);
00738         return write_float_attr(attr_name, value);
00739 }
00740 
00741 int HdfIO::write_float_attr(const string & attr_name, float value)
00742 {
00743         int err = -1;
00744         delete_attr(attr_name);
00745         hid_t dataspace = H5Screate(H5S_SCALAR);
00746         hid_t attr =
00747                 H5Acreate(cur_dataset, attr_name.c_str(), H5T_NATIVE_FLOAT, dataspace, H5P_DEFAULT);
00748 
00749         if (attr >= 0) {
00750                 err = H5Awrite(attr, H5T_NATIVE_FLOAT, &value);
00751         }
00752 
00753         H5Aclose(attr);
00754         H5Sclose(dataspace);
00755 
00756         if (err < 0) {
00757                 return 1;
00758         }
00759 
00760         return 0;
00761 }
00762 
00763 int HdfIO::write_string_attr(int image_index, const string & attr_name,
00764                                                          const string & value)
00765 {
00766         set_dataset(image_index);
00767 
00768         int err = -1;
00769         delete_attr(attr_name);
00770 
00771         hid_t datatype = H5Tcopy(H5T_C_S1);
00772         H5Tset_size(datatype, value.size() + 1);
00773         hid_t dataspace = H5Screate(H5S_SCALAR);
00774 
00775         hid_t attr = H5Acreate(cur_dataset, attr_name.c_str(), datatype, dataspace, H5P_DEFAULT);
00776         if (attr >= 0) {
00777                 err = H5Awrite(attr, datatype, (void *) value.c_str());
00778         }
00779 
00780         H5Aclose(attr);
00781         H5Sclose(dataspace);
00782         H5Tclose(datatype);
00783 
00784         if (err < 0) {
00785                 return 1;
00786         }
00787 
00788         return 0;
00789 }
00790 
00791 int HdfIO::write_array_attr(int image_index, const string & attr_name,
00792                                                         int nitems, void *data, DataType type)
00793 {
00794         if (nitems <= 0) {
00795                 return 1;
00796         }
00797         if (!data) {
00798                 throw NullPointerException("array data is NULL");
00799         }
00800 
00801         set_dataset(image_index);
00802         delete_attr(attr_name);
00803 
00804         if (type != INT && type != FLOAT) {
00805                 fprintf(stderr, "can only write INTEGER and FLOAT array");
00806                 return 1;
00807         }
00808 
00809         int err = -1;
00810         hsize_t dim[1];
00811         int perm[1];
00812         dim[0] = nitems;
00813         hsize_t sdim[] = { 1 };
00814 
00815         hid_t datatype = -1;
00816 
00817         if (type == INT) {
00818                 datatype = H5Tarray_create(H5T_NATIVE_INT, 1, dim, perm);
00819         }
00820         else if (type == FLOAT) {
00821                 datatype = H5Tarray_create(H5T_NATIVE_FLOAT, 1, dim, perm);
00822         }
00823 
00824         hid_t dataspace = H5Screate_simple(1, sdim, NULL);
00825         hid_t attr = H5Acreate(cur_dataset, attr_name.c_str(), datatype, dataspace, H5P_DEFAULT);
00826         if (attr >= 0) {
00827                 err = H5Awrite(attr, datatype, data);
00828         }
00829 
00830         H5Tclose(datatype);
00831         H5Sclose(dataspace);
00832         H5Aclose(attr);
00833 
00834         if (err < 0) {
00835                 return 1;
00836         }
00837 
00838         return 0;
00839 }
00840 
00841 
00842 int HdfIO::write_global_int_attr(const string & attr_name, int value)
00843 {
00844         hid_t tmp_dataset = cur_dataset;
00845         cur_dataset = group;
00846         int err = write_int_attr(attr_name, value);
00847         cur_dataset = tmp_dataset;
00848         return err;
00849 }
00850 
00851 #if 0
00852 int HdfIO::write_euler_attr(int image_index, const string & attr_name, float value)
00853 {
00854     string full_attr_name = string(EULER_MAGIC) + attr_name;
00855     return write_float_attr(image_index, full_attr_name, value);
00856 }
00857 #endif
00858 
00859 int HdfIO::write_mapinfo_attr(int image_index, const string & attr_name, int value)
00860 {
00861         set_dataset(image_index);
00862         delete_attr(attr_name);
00863 
00864         hsize_t dim[] = { 1 };
00865         hid_t dataspace = H5Screate_simple(1, dim, NULL);
00866         hid_t attr = H5Acreate(cur_dataset, attr_name.c_str(), mapinfo_type, dataspace, H5P_DEFAULT);
00867         H5Awrite(attr, mapinfo_type, &value);
00868         H5Aclose(attr);
00869         H5Sclose(dataspace);
00870         return 0;
00871 }
00872 
00873 
00874 int HdfIO::delete_attr(int image_index, const string & attr_name)
00875 {
00876         set_dataset(image_index);
00877 
00878         hdf_err_off();
00879         int err = H5Adelete(cur_dataset, attr_name.c_str());
00880         hdf_err_on();
00881 
00882         if (err >= 0) {
00883                 return 0;
00884         }
00885         else {
00886                 return 1;
00887         }
00888 }
00889 
00890 int HdfIO::delete_attr(const string & attr_name)
00891 {
00892         hdf_err_off();
00893         int err = H5Adelete(cur_dataset, attr_name.c_str());
00894         hdf_err_on();
00895 
00896         if (err >= 0) {
00897                 return 0;
00898         }
00899         else {
00900                 return 1;
00901         }
00902 }
00903 
00904 string HdfIO::get_compound_name(int id, const string & name)
00905 {
00906         string magic = get_item_name(COMPOUND_DATA_MAGIC);
00907         char id_str[32];
00908         sprintf(id_str, "%d", id);
00909         string compound_name = magic + "." + id_str + "." + name;
00910         return compound_name;
00911 }
00912 
00913 
00914 int HdfIO::create_compound_attr(int image_index, const string & attr_name)
00915 {
00916         string cur_dataset_name = get_compound_name(image_index, attr_name);
00917         cur_image_index = -1;
00918 
00919         hsize_t dims[1];
00920         dims[0] = 1;
00921         hid_t datatype = H5Tcopy(H5T_NATIVE_INT);
00922         hid_t dataspace = H5Screate_simple(1, dims, NULL);
00923 
00924         close_cur_dataset();
00925         cur_dataset = H5Dcreate(file, cur_dataset_name.c_str(), datatype, dataspace, H5P_DEFAULT);
00926 
00927         H5Tclose(datatype);
00928         H5Sclose(dataspace);
00929         return 0;
00930 }
00931 
00932 int HdfIO::read_ctf(Ctf & ctf, int image_index)
00933 {
00934     Dict ctf_dict;
00935     int err = read_compound_dict(CTFIT, ctf_dict, image_index);
00936     if (!err) {
00937         ctf.from_dict(ctf_dict);
00938     }
00939 
00940         return err;
00941 }
00942 
00943 void HdfIO::write_ctf(const Ctf & ctf, int image_index)
00944 {
00945     Dict ctf_dict = ctf.to_dict();
00946     write_compound_dict(CTFIT, ctf_dict, image_index);
00947 }
00948 
00949 int HdfIO::read_euler_angles(Dict & euler_angles, int image_index)
00950 {
00951     int err = read_compound_dict(EULER, euler_angles, image_index);
00952     return err;
00953 }
00954 
00955 
00956 void HdfIO::write_euler_angles(const Dict & euler_angles, int image_index)
00957 {
00958     write_compound_dict(EULER, euler_angles, image_index);
00959 }
00960 
00961 int HdfIO::read_compound_dict(Nametype compound_type,
00962                               Dict & values, int image_index)
00963 {
00964     ENTERFUNC;
00965         init();
00966 
00967         int err = 0;
00968 
00969         hid_t cur_dataset_orig = cur_dataset;
00970         string cur_dataset_name = get_compound_name(image_index, get_item_name(compound_type));
00971 
00972     hdf_err_off();
00973     cur_dataset = H5Dopen(file, cur_dataset_name.c_str());
00974     hdf_err_on();
00975 
00976         if (cur_dataset < 0) {
00977                 err = 1;
00978         }
00979         else {
00980                 err = H5Aiterate(cur_dataset, 0, attr_info, &values);
00981                 if (err < 0) {
00982                         err = 1;
00983                 }
00984         }
00985 
00986         H5Dclose(cur_dataset);
00987     cur_dataset = cur_dataset_orig;
00988 
00989         cur_image_index = -1;
00990         EXITFUNC;
00991         return err;
00992 }
00993 
00994 
00995 void HdfIO::write_compound_dict(Nametype compound_type,
00996                                 const Dict & values, int image_index)
00997 {
00998         ENTERFUNC;
00999         init();
01000 
01001         //set_dataset(image_index);
01002         hid_t cur_dataset_orig = cur_dataset;
01003 
01004         string attr_name = get_item_name(compound_type);
01005         string cur_dataset_name = get_compound_name(image_index, attr_name);
01006 
01007     hdf_err_off();
01008         cur_dataset = H5Dopen(file, cur_dataset_name.c_str());
01009     hdf_err_on();
01010 
01011         if (cur_dataset < 0) {
01012                 create_compound_attr(image_index, attr_name);
01013         }
01014     else {
01015         // remove all existing attributes first
01016         Dict attr_dict;
01017         H5Aiterate(cur_dataset, 0, attr_info, &attr_dict);
01018         vector <string> attr_keys = attr_dict.keys();
01019         for (size_t i = 0; i < attr_keys.size(); i++) {
01020             H5Adelete(cur_dataset, attr_keys[i].c_str());
01021         }
01022     }
01023 
01024     // writing new attributes
01025 
01026         vector < string > keys = values.keys();
01027         for (size_t i = 0; i < keys.size(); i++) {
01028                 float v = values[keys[i]];
01029                 write_float_attr(keys[i].c_str(), v);
01030         }
01031 
01032         H5Dclose(cur_dataset);
01033         cur_dataset = cur_dataset_orig;
01034 
01035         cur_image_index = -1;
01036         EXITFUNC;
01037 }
01038 
01039 
01040 bool HdfIO::is_complex_mode()
01041 {
01042         return false;
01043 }
01044 
01045 // always big endian
01046 bool HdfIO::is_image_big_endian()
01047 {
01048         return true;
01049 }
01050 
01051 
01052 string HdfIO::get_item_name(Nametype type)
01053 {
01054         switch (type) {
01055         case ROOT_GROUP:
01056                 return "/";
01057         case CTFIT:
01058                 return "ctfit";
01059         case NUMDATASET:
01060                 return "num_dataset";
01061         case COMPOUND_DATA_MAGIC:
01062                 return "compound";
01063     case EULER:
01064         return "euler_angles";
01065         }
01066 
01067 
01068         return "unknown";
01069 }
01070 
01071 void HdfIO::hdf_err_off()
01072 {
01073         H5Eget_auto(&old_func, &old_client_data);
01074         H5Eset_auto(0, 0);
01075 }
01076 
01077 void HdfIO::hdf_err_on()
01078 {
01079         H5Eset_auto(old_func, old_client_data);
01080 }
01081 
01082 void HdfIO::close_cur_dataset()
01083 {
01084         hdf_err_off();
01085         if (cur_dataset >= 0) {
01086                 H5Dclose(cur_dataset);
01087                 cur_dataset = -1;
01088                 cur_image_index = -1;
01089         }
01090         hdf_err_on();
01091 }
01092 
01093 void HdfIO::create_enum_types()
01094 {
01095         static int enum_types_created = 0;
01096 
01097         if (!enum_types_created) {
01098 #if 0
01099 //        Transform3D::EulerType e;
01100 //              euler_type = H5Tcreate(H5T_ENUM, sizeof(Transform3D::EulerType));
01101 //
01102 //              H5Tenum_insert(euler_type, "EMAN", (e = Transform3D::EMAN, &e));
01103 //              H5Tenum_insert(euler_type, "IMAGIC", (e = Transform3D::IMAGIC, &e));
01104 //              H5Tenum_insert(euler_type, "SPIN", (e = Transform3D::SPIN, &e));
01105 //              H5Tenum_insert(euler_type, "QUATERNION", (e = Transform3D::QUATERNION, &e));
01106 //              H5Tenum_insert(euler_type, "SGIROT", (e = Transform3D::SGIROT, &e));
01107 //              H5Tenum_insert(euler_type, "SPIDER", (e = Transform3D::SPIDER, &e));
01108 //              H5Tenum_insert(euler_type, "MRC", (e = Transform3D::MRC, &e));
01109 //
01110 //              MapInfoType m;
01111 //              mapinfo_type = H5Tcreate(H5T_ENUM, sizeof(MapInfoType));
01112 //
01113 //              H5Tenum_insert(mapinfo_type, "NORMAL", (m = NORMAL, &m));
01114 //              H5Tenum_insert(mapinfo_type, "ICOS2F_FIRST_OCTANT", (m = ICOS2F_FIRST_OCTANT, &m));
01115 //              H5Tenum_insert(mapinfo_type, "ICOS2F_FULL", (m = ICOS2F_FULL, &m));
01116 //              H5Tenum_insert(mapinfo_type, "ICOS2F_HALF", (m = ICOS2F_HALF, &m));
01117 //              H5Tenum_insert(mapinfo_type, "ICOS3F_HALF", (m = ICOS3F_HALF, &m));
01118 //              H5Tenum_insert(mapinfo_type, "ICOS3F_FULL", (m = ICOS3F_FULL, &m));
01119 //              H5Tenum_insert(mapinfo_type, "ICOS5F_HALF", (m = ICOS5F_HALF, &m));
01120 //              H5Tenum_insert(mapinfo_type, "ICOS5F_FULL", (m = ICOS5F_FULL, &m));
01121 #endif
01122                 enum_types_created = 1;
01123         }
01124 }
01125 
01126 vector < int >HdfIO::get_image_indices()
01127 {
01128         return image_indices;
01129 }
01130 
01131 int HdfIO::get_hdf_dims(int image_index, int *p_nx, int *p_ny, int *p_nz)
01132 {
01133         *p_nx = read_int_attr(image_index, "nx");
01134         *p_ny = read_int_attr(image_index, "ny");
01135         *p_nz = read_int_attr(image_index, "nz");
01136 
01137         if (*p_nx == 0 || *p_ny == 0 || *p_nz == 0) {
01138                 int ndim = 0;
01139                 int *dims = read_dims(image_index, &ndim);
01140 
01141                 if (ndim != 2 && ndim != 3) {
01142                         LOGERR("only handle 2D/3D HDF5. Your file is %dD.", ndim);
01143                         if( dims )
01144                         {
01145                                 delete [] dims;
01146                                 dims = 0;
01147                         }
01148                         return 1;
01149                 }
01150                 else {
01151                         *p_nx = dims[0];
01152                         *p_ny = dims[1];
01153                         if (ndim == 3) {
01154                                 *p_nz = dims[2];
01155                         }
01156                         else {
01157                                 *p_nz = 1;
01158                         }
01159                 }
01160                 if( dims )
01161                 {
01162                         delete [] dims;
01163                         dims = 0;
01164                 }
01165         }
01166         return 0;
01167 }
01168 
01169 herr_t HdfIO::file_info(hid_t, const char *name, void *opdata)
01170 {
01171 //      loc_id = loc_id;
01172         vector < int >*image_indices = static_cast < vector < int >*>(opdata);
01173         string magic = HdfIO::get_item_name(HdfIO::COMPOUND_DATA_MAGIC);
01174         const char *magic_str = magic.c_str();
01175 
01176         if (strncmp(name, magic_str, strlen(magic_str)) != 0) {
01177                 int id = atoi(name);
01178                 image_indices->push_back(id);
01179         }
01180 
01181         return 0;
01182 }
01183 
01184 void HdfIO::create_cur_dataset(int image_index, int nx, int ny, int nz)
01185 {
01186         int ndim = 0;
01187         int dims[3];
01188 
01189         if (nz == 1) {
01190                 ndim = 2;
01191         }
01192         else {
01193                 ndim = 3;
01194         }
01195         dims[0] = nx;
01196         dims[1] = ny;
01197         dims[2] = nz;
01198 
01199         char tmp_dataset_name[32];
01200         sprintf(tmp_dataset_name, "%d", image_index);
01201         cur_image_index = image_index;
01202 
01203         hdf_err_off();
01204         cur_dataset = H5Dopen(file, tmp_dataset_name);
01205         hdf_err_on();
01206 
01207 
01208         if (cur_dataset >= 0) {
01209                 int ndim1 = 0;
01210                 int *dims1 = read_dims(image_index, &ndim1);
01211                 Assert(ndim == ndim1);
01212 
01213                 for (int i = 0; i < ndim; i++) {
01214                         Assert(dims[i] == dims1[i]);
01215                 }
01216                 if( dims1 )
01217                 {
01218                         delete [] dims1;
01219                         dims1 = 0;
01220                 }
01221         }
01222         else {
01223                 hsize_t *sdims = new hsize_t[ndim];
01224                 for (int i = 0; i < ndim; i++) {
01225                         sdims[i] = dims[i];
01226                 }
01227 
01228                 hid_t datatype = H5Tcopy(H5T_NATIVE_FLOAT);
01229                 hid_t dataspace = H5Screate_simple(ndim, sdims, NULL);
01230 
01231                 cur_dataset = H5Dcreate(file, tmp_dataset_name, datatype, dataspace, H5P_DEFAULT);
01232 
01233                 H5Tclose(datatype);
01234                 H5Sclose(dataspace);
01235 
01236                 if( sdims )
01237                 {
01238                         delete[]sdims;
01239                         sdims = 0;
01240                 }
01241 
01242                 if (cur_dataset < 0) {
01243                         throw ImageWriteException(filename, "create dataset failed");
01244                 }
01245                 else {
01246                         increase_num_dataset();
01247                         image_indices.push_back(image_index);
01248                 }
01249         }
01250 }
01251 
01252 int HdfIO::create_region_space(hid_t * p_dataspace_id, hid_t * p_memspace_id,
01253                                                            const Region * area, int nx, int ny, int nz,
01254                                                            int image_index)
01255 {
01256         Assert(p_dataspace_id);
01257         Assert(p_memspace_id);
01258         Assert(area);
01259 
01260         if (!p_dataspace_id || !p_memspace_id || !area) {
01261                 return -1;
01262         }
01263 
01264 #if H5_VERS_MINOR > 6 || (H5_VERS_MINOR == 6 && H5_VERS_RELEASE >= 4)
01265         hsize_t offset[3];
01266 #else
01267         hssize_t offset[3];
01268 #endif
01269 
01270         hsize_t count[3];
01271 
01272         int x0 = 0, y0 = 0, z0 = 0;
01273         int xlen = 0, ylen = 0, zlen = 0;
01274 
01275         EMUtil::get_region_origins(area, &x0, &y0, &z0, nz, image_index);
01276         EMUtil::get_region_dims(area, nx, &xlen, ny, &ylen, nz, &zlen);
01277 
01278         offset[0] = static_cast < hssize_t > (x0);
01279         offset[1] = static_cast < hssize_t > (y0);
01280         offset[2] = static_cast < hssize_t > (z0);
01281 
01282         count[0] = static_cast < hsize_t > (xlen);
01283         count[1] = static_cast < hsize_t > (ylen);
01284         count[2] = static_cast < hsize_t > (zlen);
01285 
01286         *p_dataspace_id = H5Dget_space(cur_dataset);
01287 
01288         int err = H5Sselect_hyperslab(*p_dataspace_id, H5S_SELECT_SET,
01289                                                                   offset, NULL, count, NULL);
01290         if (err >= 0) {
01291                 *p_memspace_id = H5Screate_simple(3, count, NULL);
01292 
01293 #if H5_VERS_MINOR > 6 || (H5_VERS_MINOR == 6 && H5_VERS_RELEASE >= 4)
01294                 hsize_t offset_out[3];
01295 #else
01296                 hssize_t offset_out[3];
01297 #endif
01298 
01299                 offset_out[0] = 0;
01300                 offset_out[1] = 0;
01301                 offset_out[2] = 0;
01302 
01303                 err = H5Sselect_hyperslab(*p_memspace_id, H5S_SELECT_SET,
01304                                                                   offset_out, NULL, count, NULL);
01305                 if (err >= 0) {
01306                         err = 0;
01307                 }
01308         }
01309 
01310         return err;
01311 }
01312 
01313 int HdfIO::get_num_dataset()
01314 {
01315     hdf_err_off();
01316     int n = read_global_int_attr(get_item_name(NUMDATASET));
01317     hdf_err_on();
01318     return n;
01319 }
01320 
01321 
01322 #endif  //EM_HDF5

Generated on Tue Jun 11 13:40:38 2013 for EMAN2 by  doxygen 1.3.9.1