00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
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
00084
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
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);
00119 H5Fclose(file);
00120 }
00121
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
00189
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
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
00366 else if(*iter == "micrograph_id") {
00367 write_string_attr(image_index, "micrograph_id", (const char *) dict["micrograph_id"]);
00368 }
00369
00370 else {
00371 EMObject attr_val = dict[*iter];
00372
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
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
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
00431
00432
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
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
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
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
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
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
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
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
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
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
01100
01101
01102
01103
01104
01105
01106
01107
01108
01109
01110
01111
01112
01113
01114
01115
01116
01117
01118
01119
01120
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
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