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