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