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