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