EMAN2
hdfio2.cpp
Go to the documentation of this file.
1/*
2 * Author: Steven Ludtke, 04/10/2003 (sludtke@bcm.edu)
3 * Copyright (c) 2000-2006 Baylor College of Medicine
4 *
5 * This software is issued under a joint BSD/GNU license. You may use the
6 * source code in this file under either license. However, note that the
7 * complete EMAN2 and SPARX software packages have some GPL dependencies,
8 * so you are responsible for compliance with the licenses of these packages
9 * if you opt to use BSD licensing. The warranty disclaimer below holds
10 * in either instance.
11 *
12 * This complete copyright notice must be included in any revised version of the
13 * source code. Additional authorship citations may be added, but existing
14 * author citations must be preserved.
15 *
16 * This program is free software; you can redistribute it and/or modify
17 * it under the terms of the GNU General Public License as published by
18 * the Free Software Foundation; either version 2 of the License, or
19 * (at your option) any later version.
20 *
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * You should have received a copy of the GNU General Public License
27 * along with this program; if not, write to the Free Software
28 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
29 *
30 * */
31
32#ifdef USE_HDF5
33
34// #define DEBUGHDF 1
35
36#include "hdfio2.h"
37#include "geometry.h"
38#include "ctf.h"
39#include "emassert.h"
40#include "transform.h"
41#include "ctf.h"
42
43#include <iostream>
44#include <cstring>
45#include <inttypes.h>
46
47#ifndef WIN32
48 #include <sys/param.h>
49#else
50 #define MAXPATHLEN (MAX_PATH * 4)
51#endif //WIN32
52
53// Some bugs with using stdint.h, so defining our own limits. Handling as float to avoid some math mishandling
54const float INT8_min = -128.0f;
55const float INT16_min = -32768.0f;
56//const int INT32_min = -2147483648; // some potential issues using this, but we aren't supporting them at the moment anyway
57
58const float INT8_max = 127.0f;
59const float INT16_max = 32767.0f;
60//const int INT32_max = 2147483647;
61
62const float UINT8_max = 255.0f;
63const float UINT16_max = 65535.0f;
64//const unsigned int UINT32_max = 4294967295U;
65
66using namespace EMAN;
67
68static const int ATTR_NAME_LEN = 128;
69
70HdfIO2::HdfIO2(const string & fname, IOMode rw)
71: ImageIO(fname, rw), nx(1), ny(1), nz(1), is_exist(false),
72 file(-1), group(-1),
73 rendermin(0.0), rendermax(0.0), renderbits(16), renderlevel(1)
74{
75 H5dont_atexit();
76 accprop=H5Pcreate(H5P_FILE_ACCESS);
77
78 //STDIO file driver has 2G size limit on 32 bit Linux system
79 H5Pset_fapl_sec2( accprop );
80 // 0.75 ->H5D_CHUNK_CACHE_W0_DEFAULT but raises an error
81 H5Pset_cache(accprop, 0, 256, 1024*4096, 0.75); // meaningless for non-chunked data, sets the default chunk cache size per data set to 4 MB
82 //H5Pset_fapl_stdio( accprop );
83
84 // This says it's ok to use a version of the file compatible only with versions 1.8+ of the library
85 // was supposed to improve speed on attribute read/write, but made it 5x slower !?
86 //H5Pset_libver_bounds(accprop, H5F_LIBVER_18, H5F_LIBVER_LATEST);
87
88// H5Pset_fapl_core( accprop, 1048576, 0 );
89
90 hsize_t dims=1;
91 simple_space=H5Screate_simple(1,&dims,NULL);
92
93 meta_attr_dict = Dict();
94}
95
96HdfIO2::~HdfIO2()
97{
98 H5Sclose(simple_space);
99 H5Pclose(accprop);
100 if (group >= 0) {
101 H5Gclose(group);
102 }
103
104 if (file >= 0) {
105 H5Fflush(file,H5F_SCOPE_GLOBAL); // If there were no resource leaks, this wouldn't be necessary...
106 H5Fclose(file);
107 }
108
109#ifdef DEBUGHDF
110 printf("HDF: close\n");
111#endif
112}
113
114// This reads an already opened attribute and returns the results as an EMObject
115// The attribute is not closed
116 EMObject HdfIO2::read_attr(hid_t attr) {
117 hid_t type = H5Aget_type(attr);
118 hid_t spc = H5Aget_space(attr);
119 H5T_class_t cls = H5Tget_class(type);
120 size_t sz = H5Tget_size(type); // storage size, arrays handled in the 'space'
121 hssize_t pts = H5Sget_simple_extent_npoints(spc); // number of points > 1 if an array of floats or integers
122
123 EMObject ret(0);
124 char c;
125 int i;
126 float f,*fa;
127 int * ia;
128 double d;
129 char *s;
130 vector <float> fv((size_t)pts);
131 vector <int> iv((size_t)pts);
132
133 float *matrix;
134 Transform* t;
135 Ctf* ctf;
136
137 switch (cls) {
138 case H5T_INTEGER:
139 if (sz == 1) {
140 H5Aread(attr,H5T_NATIVE_CHAR,&c);
141 ret = EMObject(bool((c=='T')));
142 }
143 else if (sz == 4) {
144 if (pts == 1) {
145 H5Aread(attr,H5T_NATIVE_INT,&i);
146 ret=EMObject(i);
147 }
148 else {
149 ia=(int *)malloc((size_t)pts*sizeof(int));
150 H5Aread(attr,H5T_NATIVE_INT,ia);
151
152 for (i=0; i<pts; i++) iv[i] = ia[i];
153
154 free(ia);
155
156 ret = EMObject(iv);
157 }
158 }
159
160 break;
161 case H5T_FLOAT:
162 if (sz == 4) {
163 if (pts == 1) {
164 H5Aread(attr,H5T_NATIVE_FLOAT,&f);
165
166 ret=EMObject(f);
167 }
168 else {
169 fa = (float *)malloc((size_t)pts*sizeof(float));
170 H5Aread(attr,H5T_NATIVE_FLOAT,fa);
171
172 for (i=0; i<pts; i++) fv[i] = fa[i];
173
174 free(fa);
175
176 ret=EMObject(fv);
177 }
178 }
179 else if (sz == 8) {
180 H5Aread(attr,H5T_NATIVE_DOUBLE,&d);
181
182 ret = EMObject(d);
183 }
184
185 break;
186 case H5T_STRING:
187 s = (char *)malloc(sz+1);
188
189 H5Aread(attr,type,s);
190// H5Aread(attr,H5T_NATIVE_CHAR,s);
191
192 // It appears that CTF objects are stored in EMObjects as strings, not object pointers. Previously there was a double
193 // conversion happening, costing a LOT of speed.
194// if (s[0] == 'O' && isdigit(s[1])) {
195// ctf = new EMAN1Ctf();
196// try {
197// ctf->from_string(string(s));
198// ret = EMObject(ctf);
199// }
200// catch(...) {
201// ret=EMObject(s);
202// }
203//
204// delete ctf;
205// }
206// else if (s[0] == 'E' && isdigit(s[1])) {
207// ctf = new EMAN2Ctf();
208// try {
209// ctf->from_string(string(s));
210// ret = EMObject(ctf);
211// }
212// catch(...) {
213// ret=EMObject(s);
214// }
215//
216// delete ctf;
217// }
218// else {
219 ret = EMObject(s);
220 if ((s[0] == 'O' && isdigit(s[1])) || (s[0] == 'E' && isdigit(s[1]))) ret.force_CTF();
221// }
222
223 free(s);
224
225 break;
226 case H5T_COMPOUND:
227 matrix = (float*)malloc(12*sizeof(float));
228 H5Aread(attr, type, matrix);
229
230 t = new Transform(matrix);
231 ret = EMObject(t);
232 free(matrix);
233 delete t; t=0;
234
235 break;
236 default:
237 LOGERR("Unhandled HDF5 metadata %d", cls);
238 }
239
240 H5Sclose(spc);
241 H5Tclose(type);
242
243 return ret;
244}
245
246// This writes an attribute with specified name to a given open object
247// The attribute is opened and closed. returns 0 on success
248int HdfIO2::write_attr(hid_t loc,const char *name,EMObject obj) {
249 hid_t type=0;
250 hid_t spc=0;
251 hsize_t dims=1;
252 vector <float> fv;
253 vector <int> iv;
254
255 switch(obj.get_type())
256 {
257 case EMObject::BOOL:
258 type=H5Tcopy(H5T_NATIVE_CHAR);
259 spc=H5Scopy(simple_space);
260 break;
261 case EMObject::SHORT:
262 case EMObject::INT:
263 type=H5Tcopy(H5T_NATIVE_INT);
264 spc=H5Scopy(simple_space);
265 break;
266 case EMObject::UNSIGNEDINT:
267 type=H5Tcopy(H5T_NATIVE_UINT);
268 spc=H5Scopy(simple_space);
269 break;
270 case EMObject::FLOAT:
271 type=H5Tcopy(H5T_NATIVE_FLOAT);
272 spc=H5Scopy(simple_space);
273 break;
274 case EMObject::DOUBLE:
275 type=H5Tcopy(H5T_NATIVE_DOUBLE);
276 spc=H5Scopy(simple_space);
277 break;
278 case EMObject::STRING:
279 case EMObject::CTF:
280 type=H5Tcopy(H5T_C_S1);
281 H5Tset_size(type,strlen((const char *)obj)+1);
282 spc=H5Screate(H5S_SCALAR);
283 break;
284 case EMObject::FLOATARRAY:
285 type=H5Tcopy(H5T_NATIVE_FLOAT);
286 fv=obj;
287 dims=fv.size();
288 spc=H5Screate_simple(1,&dims,NULL);
289 break;
290 case EMObject::INTARRAY:
291 type=H5Tcopy(H5T_NATIVE_INT);
292 iv=obj;
293 dims=iv.size();
294 spc=H5Screate_simple(1,&dims,NULL);
295 break;
296 case EMObject::TRANSFORM:
297 type = H5Tcreate(H5T_COMPOUND, 12 * sizeof(float)); // Transform is a 3x4 matrix
298 H5Tinsert(type, "00", 0, H5T_NATIVE_FLOAT);
299 H5Tinsert(type, "01", 1*sizeof(float), H5T_NATIVE_FLOAT);
300 H5Tinsert(type, "02", 2*sizeof(float), H5T_NATIVE_FLOAT);
301 H5Tinsert(type, "03", 3*sizeof(float), H5T_NATIVE_FLOAT);
302 H5Tinsert(type, "10", 4*sizeof(float), H5T_NATIVE_FLOAT);
303 H5Tinsert(type, "11", 5*sizeof(float), H5T_NATIVE_FLOAT);
304 H5Tinsert(type, "12", 6*sizeof(float), H5T_NATIVE_FLOAT);
305 H5Tinsert(type, "13", 7*sizeof(float), H5T_NATIVE_FLOAT);
306 H5Tinsert(type, "20", 8*sizeof(float), H5T_NATIVE_FLOAT);
307 H5Tinsert(type, "21", 9*sizeof(float), H5T_NATIVE_FLOAT);
308 H5Tinsert(type, "22", 10*sizeof(float), H5T_NATIVE_FLOAT);
309 H5Tinsert(type, "23", 11*sizeof(float), H5T_NATIVE_FLOAT);
310 H5Tpack(type);
311
312 dims = 1; // one compound type
313 spc = H5Screate_simple(1, &dims, NULL);
314 break;
315 case EMObject::TRANSFORMARRAY:
316 case EMObject::STRINGARRAY:
317 case EMObject::EMDATA:
318 case EMObject::XYDATA:
319 case EMObject::FLOAT_POINTER:
320 case EMObject::INT_POINTER:
321 case EMObject::VOID_POINTER:
322 return -1;
323 break;
324 case EMObject::UNKNOWN:
325 break;
326 }
327
328 // we need this delete attribute call here, even we called erase_header()
329 // at the beginning of write_header(), since the "imageid_max" need be updated correctly.
330 if (H5Adelete(loc,name) < 0) {
331#ifdef DEBUGHDF
332 LOGERR("HDF: Attribute %s deletion error in write_attr().\n", name);
333#endif
334 }
335 else {
336#ifdef DEBUGHDF
337 printf("HDF: delete attribute %s successfully in write_attr().\n", name);
338#endif
339 }
340
341 hid_t attr = H5Acreate(loc,name,type,spc,H5P_DEFAULT);
342
343 bool b;
344 char c;
345 int i;
346 short si;
347 float f,*fa;
348 int * ia;
349 unsigned int ui;
350 double d;
351 const char *s;
352 Transform * tp;
353
354 switch(obj.get_type()) {
355 case EMObject::BOOL:
356 b = (bool)obj;
357
358 c = (b ? 'T' : 'F');
359
360 H5Awrite(attr,type,&c);
361 break;
362 case EMObject::SHORT:
363 si = (short)obj;
364 i = (int)si;
365 H5Awrite(attr,type,&i);
366 break;
367 case EMObject::INT:
368 i=(int)obj;
369 H5Awrite(attr,type,&i);
370 break;
371 case EMObject::UNSIGNEDINT:
372 ui=(unsigned int)obj;
373 H5Awrite(attr,type,&ui);
374 break;
375 case EMObject::FLOAT:
376 f=(float)obj;
377 H5Awrite(attr,type,&f);
378 break;
379 case EMObject::DOUBLE:
380 d=(double)obj;
381 H5Awrite(attr,type,&d);
382 break;
383 case EMObject::STRING:
384 case EMObject::CTF:
385 s=(const char *)obj;
386 H5Awrite(attr,type,s);
387 break;
388 case EMObject::FLOATARRAY:
389 fa=(float *)malloc(fv.size()*sizeof(float));
390 for (ui=0; ui<fv.size(); ui++) fa[ui]=fv[ui];
391 H5Awrite(attr,type,fa);
392 free(fa);
393 break;
394 case EMObject::INTARRAY:
395 ia=(int *)malloc(iv.size()*sizeof(int));
396 for (ui=0; ui<iv.size(); ui++) ia[ui]=iv[ui];
397 H5Awrite(attr,type,ia);
398 free(ia);
399 break;
400 case EMObject::TRANSFORM:
401 {
402 tp = (Transform *)obj;
403 fa = (float *)malloc(12*sizeof(float));
404 int r, c, k=0;
405
406 for (r=0; r<3; ++r) {
407 for (c=0; c<4; ++c) {
408 fa[k] = tp->at(r,c);
409 ++k;
410 }
411 }
412
413 H5Awrite(attr,type,fa);
414 free(fa);
415 }
416 break;
417// case EMObject::STRINGARRAY:
418// case EMObject::EMDATA:
419// case EMObject::XYDATA:
420// return -1;
421// break;
422 default:
423 LOGERR("Unhandled HDF5 metadata '%s'", name);
424
425 }
426
427 herr_t ret1 = H5Tclose(type);
428 herr_t ret2 = H5Sclose(spc);
429 herr_t ret3 = H5Aclose(attr);
430
431 if (ret1>=0 && ret2>=0 && ret3>=0) {
432 return 0;
433 }
434 else {
435 LOGERR("close error in write_attr()\n");
436 return -1;
437 }
438}
439
440// Initializes the file for read-only or read-write access
441// Data is stored under /MDF/images
442// An attribute named imageid_max stores the number of the highest
443// numbered image in the file.
444// A group is then made for each individual image, all metadata for the
445// individual images is currently associated with the GROUP, not the dataset
446// dataset-specific data could also be associated with the dataset in
447// future. At the moment, there is only a single dataset in each group.
448void HdfIO2::init()
449{
450 ENTERFUNC;
451
452 if (initialized) {
453 return;
454 }
455
456#ifdef DEBUGHDF
457 printf("HDF: init\n");
458#endif
459
460 H5Eset_auto(0, 0); // Turn off console error logging.
461
462 if (rw_mode == READ_ONLY) {
463 file = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, accprop);
464 if (file < 0) throw FileAccessException(filename);
465 }
466 else {
467 file = H5Fopen(filename.c_str(), H5F_ACC_RDWR, accprop);
468 if (file < 0) {
469 file = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, accprop);
470 if (file < 0) {
471 throw FileAccessException(filename);
472 }
473 else {
474#ifdef DEBUGHDF
475 printf("HDF: File truncated or new file created\n");
476#endif
477 }
478 }
479 }
480
481 group=H5Gopen(file,"/MDF/images");
482 if (group<0) {
483 if (rw_mode == READ_ONLY) throw ImageReadException(filename,"HDF5 file has no image data (no /MDF group)");
484 group=H5Gcreate(file,"/MDF",64); // create the group for Macromolecular data
485 if (group<0) throw ImageWriteException(filename,"Unable to add image group (/MDF) to HDF5 file");
486 H5Gclose(group);
487 group=H5Gcreate(file,"/MDF/images",4096); // create the group for images/volumes
488 if (group<0) throw ImageWriteException(filename,"Unable to add image group (/MDF/images) to HDF5 file");
489 write_attr(group,"imageid_max",EMObject(-1));
490 }
491
492 // FIXME - This section was added by Grant, presumably because DirectElectron was storing metadata items
493 // associated with the entire image group, so he automatically calls them "DDD".*, but this doesn't
494 // seem a good permanent solution...
495 else { // read the meta attributes for all images
496 int nattr=H5Aget_num_attrs(group);
497
498 char name[ATTR_NAME_LEN];
499 for (int i=0; i<nattr; i++) {
500// hid_t attr=H5Aopen_idx(group, i);
501 hid_t attr=H5Aopen_by_idx(group,".",H5_INDEX_CRT_ORDER,H5_ITER_NATIVE,i,H5P_DEFAULT,H5P_DEFAULT);
502 H5Aget_name(attr,127,name);
503
504 if (strcmp(name,"imageid_max")!=0) {
505 EMObject val=read_attr(attr);
506 meta_attr_dict["DDD."+string(name)]=val;
507 }
508
509 H5Aclose(attr);
510 }
511
512 }
513 initialized = true;
514 EXITFUNC;
515}
516
517// If this version of init() returns -1, then we have an old-style HDF5 file
518int HdfIO2::init_test()
519{
520 ENTERFUNC;
521
522 if (initialized) {
523 return 1;
524 }
525
526#ifdef DEBUGHDF
527 printf("HDF: init_test\n");
528#endif
529
530 H5Eset_auto(0, 0); // Turn off console error logging.
531
532 hid_t fileid = H5Fopen(filename.c_str(), H5F_ACC_RDONLY, H5Pcreate(H5P_FILE_ACCESS));
533 hid_t groupid = H5Gopen(fileid, "/");
534 hid_t attid = H5Aopen_name(groupid, "num_dataset");
535
536 if (attid < 0) {
537 H5Gclose(groupid);
538 H5Fclose(fileid);
539 init();
540 EXITFUNC;
541 return 0;
542 }
543 else {
544 H5Aclose(attid);
545 H5Gclose(groupid);
546 H5Fclose(fileid);
547 EXITFUNC;
548 return -1;
549 }
550}
551
552bool HdfIO2::is_valid(const void *first_block)
553{
554 ENTERFUNC;
555
556 if (first_block) {
557 char signature[8] = { char(137),char(72),char(68),char(70),char(13),char(10),char(26), char(10) };
558 if (strncmp((const char *)first_block,signature,8)==0) return true;
559
560 // const char* f=(const char *)first_block;
561 // printf("HDF: bad hdf signature %d %d %d %d %d %d %d %d",f[0],f[1],f[2],f[3],f[4],f[5],f[6],f[7]);
562
563 return false;
564 }
565 EXITFUNC;
566 return false;
567}
568
569herr_t h5_iter_attr_read( hid_t igrp, const char *attr_name, const H5A_info_t *ainfo, void *user) {
570 hid_t attr=H5Aopen_by_name(igrp,".",attr_name, H5P_DEFAULT, H5P_DEFAULT);
571
572 Dict *dict = (Dict *)user;
573
574 if (strncmp(attr_name,"EMAN.",5)==0) {
575
576 try {
577 EMObject val=HdfIO2::read_attr(attr);
578 (*dict)[attr_name+5]=val;
579 }
580 catch(...) {
581 printf("HDF: Error reading HDF attribute %s\n",attr_name+5);
582 }
583 }
584 H5Aclose(attr);
585 return 0;
586}
587
588herr_t h5_iter_attr_del( hid_t igrp, const char *attr_name, const H5A_info_t *ainfo, void *user) {
589 hid_t attr=H5Adelete_by_name(igrp,".",attr_name, H5P_DEFAULT);
590
591 return 0;
592}
593
594
595// Reads all of the attributes from the /MDF/images/<imgno> group
596int HdfIO2::read_header(Dict & dict, int image_index, const Region * area, bool)
597{
598 ENTERFUNC;
599 init();
600
601 /* Copy the meta attributes stored in /MDF/images */
602 size_t meta_attr_size = meta_attr_dict.size();
603
604 if (meta_attr_size!=0) {
605 for (size_t i=0; i<meta_attr_size; ++i) {
606 dict[meta_attr_dict.keys()[i]] = meta_attr_dict.values()[i];
607 }
608 }
609
610#ifdef DEBUGHDF
611 printf("HDF: read_head %d\n", image_index);
612#endif
613
614 int i;
615
616 // Each image is in a group for later expansion. Open the group
617 char ipath[50];
618 sprintf(ipath,"/MDF/images/%d", image_index);
619 hid_t igrp=H5Gopen(file, ipath);
620
621 if (igrp<0) {
622 char msg[40];
623 sprintf(msg,"Image %d does not exist",image_index); // yes, sprintf(), terrible I know
624 throw ImageReadException(filename,msg);
625 }
626
627 hsize_t n=0;
628 if (H5Aiterate2(igrp,H5_INDEX_NAME, H5_ITER_NATIVE, &n, &h5_iter_attr_read, &dict)) {
629 printf("Error on HDF5 iter attr\n");
630 }
631
632// int nattr=H5Aget_num_attrs(igrp);
633// char name[ATTR_NAME_LEN];
634//
635// for (i=0; i<nattr; i++) {
636// //hid_t attr=H5Aopen_idx(igrp, i);
637// hid_t attr=H5Aopen_by_idx(igrp,".",H5_INDEX_NAME,H5_ITER_NATIVE,i,H5P_DEFAULT,H5P_DEFAULT);
638// H5Aget_name(attr,127,name);
639//
640// if (strncmp(name,"EMAN.", 5)!=0) {
641// H5Aclose(attr);
642// continue;
643// }
644//
645// try {
646// EMObject val=read_attr(attr);
647// dict[name+5]=val;
648// }
649// catch(...) {
650// printf("HDF: Error reading HDF attribute %s\n",name+5);
651// }
652//
653// H5Aclose(attr);
654// }
655
656 if (dict.has_key("ctf")) {
657 dict["ctf"].force_CTF();
658// string ctfString = (string)dict["ctf"];
659//
660// if (ctfString.substr(0, 1) == "O") {
661// Ctf * ctf_ = new EMAN1Ctf();
662// ctf_->from_string(ctfString);
663// dict.erase("ctf");
664// dict["ctf"] = ctf_;
665// delete ctf_;
666// }
667// else if (ctfString.substr(0, 1) == "E") {
668// Ctf * ctf_ = new EMAN2Ctf();
669// ctf_->from_string(ctfString);
670// dict.erase("ctf");
671// dict["ctf"] = ctf_;
672// delete ctf_;
673// }
674 }
675
676 if (area) {
677 check_region(area, IntSize(dict["nx"], dict["ny"], dict["nz"]), false, false);
678
679 dict["nx"] = area->get_width();
680 dict["ny"] = area->get_height();
681 dict["nz"] = area->get_depth();
682
683 if (dict.has_key("apix_x") && dict.has_key("apix_y") && dict.has_key("apix_z"))
684 {
685 if (dict.has_key("origin_x") && dict.has_key("origin_y") && dict.has_key("origin_z"))
686 {
687 float xorigin = dict["origin_x"];
688 float yorigin = dict["origin_y"];
689 float zorigin = dict["origin_z"];
690
691 float apix_x = dict["apix_x"];
692 float apix_y = dict["apix_y"];
693 float apix_z = dict["apix_z"];
694
695 dict["origin_x"] = xorigin + apix_x * area->origin[0];
696 dict["origin_y"] = yorigin + apix_y * area->origin[1];
697 dict["origin_z"] = zorigin + apix_z * area->origin[2];
698 }
699 }
700 }
701
702 H5Gclose(igrp);
703
704 //Get the data type from data set, HDF5 file header attribute 'datatype' may be wrong
705 sprintf(ipath,"/MDF/images/%d/image",image_index);
706 hid_t ds=H5Dopen(file,ipath);
707
708 // FIXME - This isn't valid any more, since we support signed and unsigned, and have compression
709 if (ds > 0) { // ds > 0 means successfully opened the dataset
710 hid_t dt = H5Dget_type(ds);
711
712 switch(H5Tget_size(dt)) {
713 case 4:
714 dict["datatype"] = (int)EMUtil::EM_FLOAT;
715 break;
716 case 2:
717 dict["datatype"] = (int)EMUtil::EM_USHORT;
718 break;
719 case 1:
720 dict["datatype"] = (int)EMUtil::EM_UCHAR;
721 break;
722 default:
723 throw ImageReadException(filename, "EMAN does not support this data type.");
724 }
725
726 H5Tclose(dt);
727 }
728
729 H5Dclose(ds);
730
731 EXITFUNC;
732 return 0;
733}
734
735// This erases any existing attributes from the image group
736// prior to writing a new header. For a new image there
737// won't be any, so this should be harmless.
738int HdfIO2::erase_header(int image_index)
739{
740 ENTERFUNC;
741
742 if (image_index < 0) return 0; // image_index<0 for appending image, no need for erasing
743
744 init();
745
746#ifdef DEBUGHDF
747 printf("HDF: erase_head %d\n",image_index);
748#endif
749
750 int i;
751
752 // Each image is in a group for later expansion. Open the group
753 char ipath[50];
754 sprintf(ipath,"/MDF/images/%d", image_index);
755 hid_t igrp=H5Gopen(file, ipath);
756
757
758 hsize_t n=0;
759 if (H5Aiterate2(igrp,H5_INDEX_NAME, H5_ITER_NATIVE, &n, &h5_iter_attr_del, NULL)) {
760 printf("Error on HDF5 erase_header\n");
761 }
762
763
764// int nattr=H5Aget_num_attrs(igrp);
765//
766// char name[ATTR_NAME_LEN];
767//
768// for (i=0; i<nattr; i++) {
769// // hid_t attr = H5Aopen_idx(igrp, 0); // use 0 as index here, since the H5Adelete() will change the index
770// // H5Aget_name(attr,127,name);
771// // H5Aclose(attr);
772// //
773// // if (H5Adelete(igrp,name) < 0) {
774// // LOGERR("attribute %s deletion error in erase_header().\n", name);
775// // }
776// hid_t attr=H5Adelete_by_idx(igrp,".",H5_INDEX_UNKNOWN,H5_ITER_NATIVE,i,H5P_DEFAULT);
777//
778// }
779
780 H5Gclose(igrp);
781 EXITFUNC;
782 return 0;
783}
784
785// TODO : incomplete
786int HdfIO2::read_data_8bit(unsigned char *data, int image_index, const Region *area, bool is_3d, float minval, float maxval) {
787 ENTERFUNC;
788#ifdef DEBUGHDF
789 printf("HDF: read_data_8bit %d\n",image_index);
790#endif
791
792 char ipath[50];
793 sprintf(ipath,"/MDF/images/%d/image",image_index);
794 hid_t ds = H5Dopen(file,ipath);
795
796 if (ds < 0) throw ImageWriteException(filename,"Image does not exist");
797
798 hid_t spc=H5Dget_space(ds);
799 hid_t dt = H5Dget_type(ds);
800
801 hsize_t dims_out[3];
802 hsize_t rank = H5Sget_simple_extent_ndims(spc);
803
804 H5Sget_simple_extent_dims(spc, dims_out, NULL);
805
806 if (rank == 1) {
807 nx = dims_out[0];
808 ny = 1;
809 nz = 1;
810 }
811 else if (rank == 2) {
812 nx = dims_out[1];
813 ny = dims_out[0];
814 nz = 1;
815 }
816 else if (rank == 3) {
817 nx = dims_out[2];
818 ny = dims_out[1];
819 nz = dims_out[0];
820 }
821
822 if (area) {
823 hid_t memoryspace = 0;
824
825 /* Get the file dataspace - the region we want to read in the file */
826 int x0 = 0, y0 = 0, z0 = 0; // the coordinates for up left corner, trim to be within image bound
827 int x1 = 0, y1 = 0, z1 = 0; // the coordinates for down right corner, trim to be within image bound
828 int nx1 = 1, ny1 = 1, nz1 = 1; // dimensions of the sub-region, actual region read from file
829
830 if (rank == 3) {
831 hsize_t doffset[3]; /* hyperslab offset in the file */
832
833 doffset[2] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
834 doffset[1] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
835 doffset[0] = (hsize_t)(area->z_origin() < 0 ? 0 : area->z_origin());
836
837 x0 = (int)doffset[0];
838 y0 = (int)doffset[1];
839 z0 = (int)doffset[2];
840
841 z1 = (int)(area->x_origin() + area->get_width());
842 z1 = (int)(z1 > static_cast<int>(nx) ? nx : z1);
843
844 y1 = (int)(area->y_origin() + area->get_height());
845 y1 = (int)(y1 > static_cast<int>(ny) ? ny : y1);
846
847 if (y1 <= 0) {
848 y1 = 1;
849 }
850
851 x1 = (int)(area->z_origin() + area->get_depth());
852
853 x1 = (int)(x1 > static_cast<int>(nz) ? nz : x1);
854 if (x1 <= 0) {
855 x1 = 1;
856 }
857
858 if (x1 < x0 || y1< y0 || z1 < z0) return 0; //out of bounds, this is fine, nothing happens
859
860 hsize_t dcount[3]; /* size of the hyperslab in the file */
861
862 dcount[0] = x1 - doffset[0];
863 dcount[1] = y1 - doffset[1];
864 dcount[2] = z1 - doffset[2];
865
866 H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL,
867 (const hsize_t*)dcount, NULL);
868
869 /* Define memory dataspace - the memory we will created for the region */
870 hsize_t dims[3]; /* size of the region in the memory */
871
872 dims[0] = dcount[2]?dcount[2]:1;
873 dims[1] = dcount[1]?dcount[1]:1;
874 dims[2] = dcount[0]?dcount[0]:1;
875
876 nx1 = (int)dims[0];
877 ny1 = (int)dims[1];
878 nz1 = (int)dims[2];
879
880 memoryspace = H5Screate_simple(3, dims, NULL);
881 }
882 else if (rank == 2) {
883 hsize_t doffset[2]; /* hyperslab offset in the file */
884
885 doffset[1] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
886 doffset[0] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
887
888 x0 = (int)doffset[0];
889 y0 = (int)doffset[1];
890 z0 = 1;
891
892 y1 = (int)(area->x_origin() + area->get_width());
893 y1 = (int)(y1 > static_cast<int>(nx) ? nx : y1);
894
895 x1 = (int)(area->y_origin() + area->get_height());
896 x1 = (int)(x1 > static_cast<int>(ny) ? ny : x1);
897
898 if (x1 <= 0) {
899 x1 = 1;
900 }
901
902 z1 = 1;
903
904 if (x1 < x0 || y1< y0) return 0; // out of bounds, this is fine, nothing happens
905
906 hsize_t dcount[2]; /* size of the hyperslab in the file */
907 dcount[0] = x1 - doffset[0];
908 dcount[1] = y1 - doffset[1];
909
910 H5Sselect_none(spc);
911 H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL,
912 (const hsize_t*)dcount, NULL);
913
914 /* Define memory dataspace - the memory we will created for the region */
915 hsize_t dims[2]; /* size of the region in the memory */
916
917 dims[0] = (hsize_t)(dcount[1]?dcount[1]:1);
918 dims[1] = (hsize_t)(dcount[0]?dcount[0]:1);
919
920 nx1 = (int)dims[0];
921 ny1 = (int)dims[1];
922 nz1 = 1;
923
924 memoryspace = H5Screate_simple(2, dims, NULL);
925 }
926
927 if ((area->x_origin()>=0) && (area->y_origin()>=0) && (area->z_origin()>=0) &&
928 ((hsize_t)(area->x_origin() + area->get_width())<=nx) &&
929 ((hsize_t)(area->y_origin() + area->get_height())<=ny) &&
930 ((hsize_t)(area->z_origin() + area->get_depth())<=nz)) { // the region is in boundary
931
932 H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,data);
933 }
934 else { // the region are partial out of boundary
935 /* When the requested region has some part out of image boundary,
936 * we need read the sub-area which is within image,
937 * and fill the out of boundary part with zero.
938 * We actually read the sub-region from HDF by hyperslab I/O,
939 * then copy it back to the pre-allocated region. */
940
941 float * subdata = new float[nx1*ny1*nz1];
942
943 H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,subdata);
944
945 int xd0=0, yd0=0, zd0=0; // The coordinates of the top-left corner sub-region in region
946 size_t clipped_row_size = 0;
947
948 if (rank == 3) {
949 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
950 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
951 zd0 = (int) (area->z_origin() < 0 ? -area->z_origin() : 0);
952 clipped_row_size = (z1-z0)* sizeof(float);
953 }
954 else if (rank == 2) {
955 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
956 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
957 clipped_row_size = (y1-y0)* sizeof(float);
958 }
959
960 int src_secsize = nx1 * ny1;
961 int dst_secsize = (int)(area->get_width())*(int)(area->get_height());
962
963 float * src = subdata;
964 unsigned char * dst = data + zd0*dst_secsize + yd0*(int)(area->get_width()) + xd0;
965
966 int src_gap = src_secsize - (y1-y0) * nx1;
967 int dst_gap = dst_secsize - (y1-y0) * (int)(area->get_width());
968
969 for (int i = 0; i<nz1; ++i) {
970 for (int j = 0; j<ny1; ++j) {
971 EMUtil::em_memcpy(dst, src, clipped_row_size);
972
973 src += nx1;
974 dst += (int)(area->get_width());
975 }
976
977 src += src_gap;
978 dst += dst_gap;
979 }
980
981 delete [] subdata;
982 }
983
984 H5Sclose(memoryspace);
985 } else {
986 hsize_t size = (hsize_t)nx*ny*nz;
987 hsize_t i=0;
988 hsize_t j=0;
989
990 unsigned short *usdata = (unsigned short *) data;
991 unsigned char *cdata = (unsigned char *) data;
992
993 switch(H5Tget_size(dt)) {
994 case 4:
995 H5Dread(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
996
997 break;
998 case 2:
999 H5Dread(ds,H5T_NATIVE_USHORT,spc,spc,H5P_DEFAULT,usdata);
1000
1001 for (i = 0; i < size; ++i) {
1002 j = size - 1 - i;
1003 data[j] = static_cast < float >(usdata[j]);
1004 }
1005
1006 break;
1007 case 1:
1008 H5Dread(ds,H5T_NATIVE_UCHAR,spc,spc,H5P_DEFAULT,cdata);
1009
1010 for (i = 0; i < size; ++i) {
1011 j = size - 1 - i;
1012 data[j] = static_cast < float >(cdata[j]);
1013 }
1014
1015 break;
1016 default:
1017 throw ImageReadException(filename, "EMAN does not support this data type.");
1018 }
1019 }
1020
1021 H5Tclose(dt);
1022 H5Sclose(spc);
1023 H5Dclose(ds);
1024
1025 EXITFUNC;
1026 return 0;
1027}
1028
1029int HdfIO2::read_data(float *data, int image_index, const Region *area, bool)
1030{
1031 ENTERFUNC;
1032#ifdef DEBUGHDF
1033 printf("HDF: read_data %d\n",image_index);
1034#endif
1035
1036 char ipath[50];
1037 sprintf(ipath,"/MDF/images/%d/image",image_index);
1038 hid_t ds = H5Dopen(file,ipath);
1039
1040 if (ds < 0) throw ImageReadException(filename,"Image does not exist");
1041
1042 hid_t spc=H5Dget_space(ds);
1043 hid_t dt = H5Dget_type(ds);
1044
1045 hsize_t dims_out[3];
1046 hsize_t rank = H5Sget_simple_extent_ndims(spc);
1047
1048 H5Sget_simple_extent_dims(spc, dims_out, NULL);
1049
1050 if (rank == 1) {
1051 nx = dims_out[0];
1052 ny = 1;
1053 nz = 1;
1054 }
1055 else if (rank == 2) {
1056 nx = dims_out[1];
1057 ny = dims_out[0];
1058 nz = 1;
1059 }
1060 else if (rank == 3) {
1061 nx = dims_out[2];
1062 ny = dims_out[1];
1063 nz = dims_out[0];
1064 }
1065 hsize_t size = (hsize_t)nx*(hsize_t)ny*(hsize_t)nz;
1066
1067 if (area) {
1068 hid_t memoryspace = 0;
1069
1070 /* Get the file dataspace - the region we want to read in the file */
1071 int x0 = 0, y0 = 0, z0 = 0; // the coordinates for up left corner, trim to be within image bound
1072 int x1 = 0, y1 = 0, z1 = 0; // the coordinates for down right corner, trim to be within image bound
1073 int nx1 = 1, ny1 = 1, nz1 = 1; // dimensions of the sub-region, actual region read from file
1074
1075 if (rank == 3) {
1076 hsize_t doffset[3]; /* hyperslab offset in the file */
1077
1078// if ((floor(area->get_width()/256.0)+2)*(floor(area->get_height()/256.0)+2)*area->get_depth()>16)
1079// printf("Region too large: %1.0f %1.0f %1.0f\n",area->get_width(),area->get_height(),area->get_depth());
1080
1081 doffset[2] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
1082 doffset[1] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
1083 doffset[0] = (hsize_t)(area->z_origin() < 0 ? 0 : area->z_origin());
1084
1085 x0 = (int)doffset[0];
1086 y0 = (int)doffset[1];
1087 z0 = (int)doffset[2];
1088
1089 z1 = (int)(area->x_origin() + area->get_width());
1090 z1 = (int)(z1 > static_cast<int>(nx) ? nx : z1);
1091
1092 y1 = (int)(area->y_origin() + area->get_height());
1093 y1 = (int)(y1 > static_cast<int>(ny) ? ny : y1);
1094
1095 if (y1 <= 0) {
1096 y1 = 1;
1097 }
1098
1099 x1 = (int)(area->z_origin() + area->get_depth());
1100
1101 x1 = (int)(x1 > static_cast<int>(nz) ? nz : x1);
1102 if (x1 <= 0) {
1103 x1 = 1;
1104 }
1105
1106 if (x1 < x0 || y1< y0 || z1 < z0) return 0; //out of bounds, this is fine, nothing happens
1107
1108 hsize_t dcount[3]; /* size of the hyperslab in the file */
1109
1110 dcount[0] = x1 - doffset[0];
1111 dcount[1] = y1 - doffset[1];
1112 dcount[2] = z1 - doffset[2];
1113
1114 H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL,
1115 (const hsize_t*)dcount, NULL);
1116
1117 /* Define memory dataspace - the memory we will created for the region */
1118 hsize_t dims[3]; /* size of the region in the memory */
1119
1120 dims[0] = dcount[2]?dcount[2]:1;
1121 dims[1] = dcount[1]?dcount[1]:1;
1122 dims[2] = dcount[0]?dcount[0]:1;
1123
1124 nx1 = (int)dims[0];
1125 ny1 = (int)dims[1];
1126 nz1 = (int)dims[2];
1127
1128 memoryspace = H5Screate_simple(3, dims, NULL);
1129 }
1130 else if (rank == 2) {
1131 hsize_t doffset[2]; /* hyperslab offset in the file */
1132
1133 doffset[1] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
1134 doffset[0] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
1135
1136 x0 = (int)doffset[0];
1137 y0 = (int)doffset[1];
1138 z0 = 1;
1139
1140 y1 = (int)(area->x_origin() + area->get_width());
1141 y1 = (int)(y1 > static_cast<int>(nx) ? nx : y1);
1142
1143 x1 = (int)(area->y_origin() + area->get_height());
1144 x1 = (int)(x1 > static_cast<int>(ny) ? ny : x1);
1145
1146 if (x1 <= 0) {
1147 x1 = 1;
1148 }
1149
1150 z1 = 1;
1151
1152 if (x1 < x0 || y1< y0) return 0; // out of bounds, this is fine, nothing happens
1153
1154 hsize_t dcount[2]; /* size of the hyperslab in the file */
1155 dcount[0] = x1 - doffset[0];
1156 dcount[1] = y1 - doffset[1];
1157
1158 H5Sselect_none(spc);
1159 H5Sselect_hyperslab (spc, H5S_SELECT_SET, (const hsize_t*)doffset, NULL,
1160 (const hsize_t*)dcount, NULL);
1161
1162 /* Define memory dataspace - the memory we will created for the region */
1163 hsize_t dims[2]; /* size of the region in the memory */
1164
1165 dims[0] = (hsize_t)(dcount[1]?dcount[1]:1);
1166 dims[1] = (hsize_t)(dcount[0]?dcount[0]:1);
1167
1168 nx1 = (int)dims[0];
1169 ny1 = (int)dims[1];
1170 nz1 = 1;
1171
1172 memoryspace = H5Screate_simple(2, dims, NULL);
1173 }
1174
1175 if ((area->x_origin()>=0) && (area->y_origin()>=0) && (area->z_origin()>=0) &&
1176 ((hsize_t)(area->x_origin() + area->get_width())<=nx) &&
1177 ((hsize_t)(area->y_origin() + area->get_height())<=ny) &&
1178 ((hsize_t)(area->z_origin() + area->get_depth())<=nz)) { // the region is in boundary
1179
1180 H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,data);
1181 }
1182 else { // the region are partial out of boundary
1183 /* When the requested region has some part out of image boundary,
1184 * we need read the sub-area which is within image,
1185 * and fill the out of boundary part with zero.
1186 * We actually read the sub-region from HDF by hyperslab I/O,
1187 * then copy it back to the pre-allocated region. */
1188
1189 float * subdata = new float[nx1*ny1*nz1];
1190
1191 H5Dread(ds,H5T_NATIVE_FLOAT,memoryspace,spc,H5P_DEFAULT,subdata);
1192
1193 int xd0=0, yd0=0, zd0=0; // The coordinates of the top-left corner sub-region in region
1194 size_t clipped_row_size = 0;
1195
1196 if (rank == 3) {
1197 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
1198 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
1199 zd0 = (int) (area->z_origin() < 0 ? -area->z_origin() : 0);
1200 clipped_row_size = (z1-z0)* sizeof(float);
1201 }
1202 else if (rank == 2) {
1203 xd0 = (int) (area->x_origin() < 0 ? -area->x_origin() : 0);
1204 yd0 = (int) (area->y_origin() < 0 ? -area->y_origin() : 0);
1205 clipped_row_size = (y1-y0)* sizeof(float);
1206 }
1207
1208 int src_secsize = nx1 * ny1;
1209 int dst_secsize = (int)(area->get_width())*(int)(area->get_height());
1210
1211 float * src = subdata;
1212 float * dst = data + zd0*dst_secsize + yd0*(int)(area->get_width()) + xd0;
1213
1214 int src_gap = src_secsize - (y1-y0) * nx1;
1215 int dst_gap = dst_secsize - (y1-y0) * (int)(area->get_width());
1216
1217 for (int i = 0; i<nz1; ++i) {
1218 for (int j = 0; j<ny1; ++j) {
1219 EMUtil::em_memcpy(dst, src, clipped_row_size);
1220
1221 src += nx1;
1222 dst += (int)(area->get_width());
1223 }
1224
1225 src += src_gap;
1226 dst += dst_gap;
1227 }
1228
1229 delete [] subdata;
1230 }
1231 size = (hsize_t) area->get_width()*(hsize_t)area->get_height()*(hsize_t)area->get_depth();
1232
1233 H5Sclose(memoryspace);
1234 } else {
1235 hsize_t size = (hsize_t)nx*ny*nz;
1236 hsize_t i=0;
1237 hsize_t j=0;
1238
1239 unsigned short *usdata = (unsigned short *) data;
1240 unsigned char *cdata = (unsigned char *) data;
1241
1242 H5Dread(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
1243
1244 }
1245
1246 H5Tclose(dt);
1247 H5Sclose(spc);
1248 H5Dclose(ds);
1249 // Rescale data on read if bit reduction took place
1250 sprintf(ipath,"/MDF/images/%d",image_index);
1251 hid_t igrp=H5Gopen(file,ipath);
1252 hid_t iattr=H5Aopen_name(igrp,"EMAN.stored_renderbits");
1253 if (iattr>=0) {
1254 renderbits=(int)read_attr(iattr);
1255// printf("stored renderbits %d\n",renderbits);
1256 H5Aclose(iattr);
1257 if (renderbits>0) {
1258 iattr=H5Aopen_name(igrp,"EMAN.stored_rendermax");
1259 if (iattr>=0) {
1260 rendermax=(float)read_attr(iattr);
1261 H5Aclose(iattr);
1262 iattr=H5Aopen_name(igrp,"EMAN.stored_rendermin");
1263 if (iattr>=0) {
1264 rendermin=(float)read_attr(iattr);
1265 H5Aclose(iattr);
1266 float RUMAX = (1<<renderbits)-1.0f;
1267 for (size_t i=0; i<size; i++) data[i]=(data[i]/RUMAX)*(rendermax-rendermin)+rendermin;
1268 }
1269 }
1270 }
1271 }
1272
1273 H5Gclose(igrp);
1274
1275 EXITFUNC;
1276 return 0;
1277}
1278
1279// Writes all attributes in 'dict' to the image group
1280// Creation of the image dataset is also handled here
1281int HdfIO2::write_header(const Dict & dict, int image_index, const Region* area,
1282 EMUtil::EMDataType dt, bool endian)
1283{
1284#ifdef DEBUGHDF
1285 printf("HDF: write_head %d\n",image_index);
1286#endif
1287
1288 ENTERFUNC;
1289
1290 init();
1291
1292 nx = (int)dict["nx"];
1293 ny = (int)dict["ny"];
1294 nz = (int)dict["nz"];
1295
1296 if (image_index<0) {
1297 image_index = get_nimg();
1298 }
1299
1300 // If image_index<0 append, and make sure the max value in the file is correct
1301 // though this is normally handled by EMData.write_image()
1302 hid_t attr=H5Aopen_name(group,"imageid_max");
1303 int nimg = read_attr(attr);
1304 H5Aclose(attr);
1305
1306 unsigned int i;
1307
1308 if (image_index < 0) image_index=nimg+1;
1309
1310 if (image_index > nimg) {
1311 write_attr(group,(const char *)"imageid_max",EMObject(image_index));
1312 }
1313
1314 // Each image is in a group for later expansion. Open the group, create if necessary
1315 char ipath[50];
1316 sprintf(ipath,"/MDF/images/%d",image_index);
1317 hid_t igrp=H5Gopen(file,ipath);
1318
1319 if (igrp < 0) { // group doesn't exist
1320 is_exist = false;
1321 // Need to create a new image group
1322 igrp=H5Gcreate(file,ipath,64); // The image is a group, with attributes on the group
1323
1324 if (igrp < 0) throw ImageWriteException(filename,"Unable to add /MDF/images/# to HDF5 file");
1325 }
1328 else {
1329 is_exist = true;
1330// int nattr=H5Aget_num_attrs(igrp);
1331// char name[ATTR_NAME_LEN];
1332 Dict dict2;
1333 hsize_t n=0;
1334 if (H5Aiterate2(igrp,H5_INDEX_NAME, H5_ITER_NATIVE, &n, &h5_iter_attr_read, &dict2)) {
1335 printf("Error on HDF5 iter attr\n");
1336 }
1337
1338// for (int i=0; i<nattr; i++) {
1339// hid_t attr=H5Aopen_idx(igrp, i);
1340// H5Aget_name(attr,127,name);
1341//
1342// if (strncmp(name,"EMAN.", 5)!=0) {
1343// H5Aclose(attr);
1344// continue;
1345// }
1346//
1347// EMObject val=read_attr(attr);
1348// dict2[name+5]=val;
1349// H5Aclose(attr);
1350//
1351// }
1352
1353 if (!dict2.has_key("datatype")) { // by default, HDF5 is written as float
1354 dict2["datatype"] = (int)EMUtil::EM_FLOAT;
1355 }
1356
1357 if (area) {
1358 check_region(area, IntSize(dict2["nx"], dict2["ny"], dict2["nz"]), false, true);
1359 }
1360 else {
1361 erase_header(image_index);
1362
1363// printf("eraseh %d %d %d %d %d %d %d %d %d\n",(int)dict["nx"],(int)dict["ny"],(int)dict["nz"],(int)dict2["nx"],(int)dict2["ny"],(int)dict2["nz"],(int)dict["datatype"],(int)dict2["datatype"],(int)dt);
1364 // change the size or data type of a image,
1365 // the existing data set is invalid, unlink it
1366 if ((int)dict["nx"]*(int)dict["ny"]*(int)dict["nz"] !=
1367 (int)dict2["nx"]*(int)dict2["ny"]*(int)dict2["nz"] ||
1368 (int)dict["datatype"] != (int)dict2["datatype"] ||
1369 dict2.has_key("render_compress_level") || (int)dt==(int)EMUtil::EM_COMPRESSED) {
1370
1371// printf("erase\n");
1372 sprintf(ipath,"/MDF/images/%d/image",image_index);
1373
1374 H5Gunlink(igrp, ipath);
1375 }
1376 }
1377 }
1378
1379 if (! area) {
1380 // Write the attributes to the group
1381 vector <string> keys=dict.keys();
1382
1383 for (i=0; i<keys.size(); i++) {
1384 // These keys are written later if used. They MUST not be copied from a previous read operation or non-compressed
1385 // images will wind up being incorrectly scaled!
1386 if (keys[i]==string("stored_rendermin") || keys[i]==string("stored_rendermax") || keys[i]==string("stored_renderbits") ||
1387 keys[i]==string("render_min") || keys[i]==string("render_max") || keys[i]==string("render_bits")) continue;
1388 string s("EMAN.");
1389 s+=keys[i];
1390 write_attr(igrp,s.c_str(),dict[keys[i]]);
1391 }
1392 }
1393
1394 H5Gclose(igrp);
1395
1396 // Set render_min and render_max from EMData attr's if possible.
1397 if (dict.has_key("render_compress_level")) renderlevel=(float)dict["render_compress_level"];
1398 EMUtil::getRenderLimits(dict, rendermin, rendermax, renderbits);
1399
1400 EXITFUNC;
1401 return 0;
1402}
1403
1404// Writes the actual image data to the corresponding dataset (already created)
1405int HdfIO2::write_data(float *data, int image_index, const Region* area,
1406 EMUtil::EMDataType dt, bool)
1407{
1408 ENTERFUNC;
1409
1410#ifdef DEBUGHDF
1411 printf("HDF: write_data %d\n",image_index);
1412#endif
1413
1414 if (image_index < 0) {
1415 hid_t attr=H5Aopen_name(group,"imageid_max");
1416 image_index = read_attr(attr);
1417 H5Aclose(attr);
1418 }
1419
1420 char ipath[50];
1421 hid_t spc; //dataspace
1422 hid_t ds; //dataset
1423
1424 sprintf(ipath, "/MDF/images/%d/image", image_index);
1425
1426 // Now create the image dataspace (not used for region writing)
1427 hsize_t rank = 0;
1428 if (nz == 1 && ny == 1) {
1429 hsize_t dims[1]= { nx };
1430 spc=H5Screate_simple(1,dims,NULL);
1431 }
1432 else if (nz == 1) {
1433 hsize_t dims[2]= { ny,nx };
1434 spc=H5Screate_simple(2,dims,NULL);
1435 }
1436 else {
1437 hsize_t dims[3]= { nz, ny, nx };
1438 spc=H5Screate_simple(3,dims,NULL);
1439 }
1440
1441 // Setup a standard HDF datatype
1442 hid_t hdt=0;
1443 switch(dt) {
1444 case EMUtil::EM_FLOAT: hdt=H5T_NATIVE_FLOAT; break;
1445 case EMUtil::EM_USHORT: hdt=H5T_NATIVE_USHORT; break;
1446 case EMUtil::EM_UCHAR: hdt=H5T_NATIVE_UCHAR; break;
1447 case EMUtil::EM_SHORT: hdt=H5T_NATIVE_SHORT; break;
1448 case EMUtil::EM_CHAR: hdt=H5T_NATIVE_CHAR; break;
1449 case EMUtil::EM_COMPRESSED:
1450 if (renderbits<=0) hdt=H5T_NATIVE_FLOAT;
1451 else if (renderbits<=8) hdt=H5T_NATIVE_UCHAR;
1452 else if (renderbits<=16) hdt=H5T_NATIVE_USHORT;
1453 else throw ImageWriteException(filename,"Bit reduced compressed HDF5 files may not use more than 16 bits. For native float, set 0 bits.");
1454 break;
1455 default:
1456 throw ImageWriteException(filename,"HDF5 does not support this data format");
1457 }
1458
1459 if (nx==1 && dt==EMUtil::EM_COMPRESSED) {
1460 printf("Warning: HDF compressed mode not supported when nx=1\n");
1461 dt=EMUtil::EM_FLOAT;
1462 }
1463
1464 // Now, we try and open the dataset
1465 ds = H5Dopen(file,ipath);
1466
1467 // In theory we will catch all of these mismatches when writing the header above, and unlinking mismatching
1468 // data objects (not truly deleting, could make files bigger). This code shouldn't cost much, so will leave
1469 // it here to double-check 11/14/20
1470 if (ds>=0) {
1471 hid_t cpl=H5Dget_create_plist(ds);
1472 H5D_layout_t layout = H5Pget_layout(cpl);
1473 H5Pclose(cpl);
1474
1475// printf("%d %d %d %d ",(int)layout,(int)(dt==EMUtil::EM_COMPRESSED),(int)H5D_CHUNKED,(int)H5D_COMPACT);
1476 if ((layout==H5D_CHUNKED && dt!=EMUtil::EM_COMPRESSED) || (layout!=H5D_CHUNKED && dt==EMUtil::EM_COMPRESSED)) {
1477 H5Sclose(spc);
1478 H5Dclose(ds);
1479 throw ImageWriteException(filename,"HDF ERROR: Trying to overwrite an image with compression mismatch");
1480 }
1481
1482 // causes issues with region writing, and not sure it's really necessary
1483// hid_t spc_file = H5Dget_space(ds);
1484// if (H5Sget_simple_extent_ndims(spc_file)!=rank) {
1485// H5Sclose(spc_file);
1486// H5Sclose(spc);
1487// H5Dclose(ds);
1488// throw ImageWriteException(filename,"HDF ERROR: Trying to overwrite an image with different dimensionality");
1489// }
1490// H5Sclose(spc_file);
1491
1492 hid_t fdt = H5Dget_type(ds);
1493// printf("%d %d\n",(int)fdt,(int)dt);
1494 if (!H5Tequal(fdt,dt)) {
1495 H5Tclose(fdt);
1496 H5Sclose(spc);
1497 H5Dclose(ds);
1498 throw ImageWriteException(filename,"HDF ERROR: Trying to overwrite an image with different data type");
1499 }
1500
1501
1502 }
1503
1504 if (ds < 0) { //new dataset
1505 hid_t plist = H5Pcreate(H5P_DATASET_CREATE); // we could just use H5P_DEFAULT for non-compressed
1506 if (dt==EMUtil::EM_COMPRESSED) {
1507 if (nz==1) {
1508 hsize_t cdims[2] = { ny>256?256:ny, nx>256?256:nx}; // slice-wise reading common in 3D so 2-D chunks
1509 H5Pset_chunk(plist,2,cdims); // uses only the first 2 elements
1510 }
1511 else {
1512 hsize_t cdims[3] = { 1, ny>256?256:ny, nx>256?256:nx}; // slice-wise reading common in 3D so 2-D chunks
1513 H5Pset_chunk(plist,3,cdims);
1514 }
1515 // H5Pset_scaleoffset(plist,H5Z_SO_FLOAT_DSCALE,2); // doesn't seem to work right?, anyway some conceptual problems
1516 // H5Pset_shuffle(plist); // rearrange bytes, seems to have zero impact, maybe deflate is internally doing something?
1517 H5Pset_deflate(plist, renderlevel); // zlib level default is 1
1518 //conclusion is that SZIP compresses roughly as well as GZIP3, but is twice as fast (for 4 bit cryoem data)
1519 // While this is good, it isn't worth the IP hassles right now. We can revisit the issue later if a good
1520 // open license library starts being widely used
1521 //int r=H5Pset_szip (plist, H5_SZIP_NN_OPTION_MASK, 32); // szip with 32 pixels per block (NN (2 stage) vs EC), NN definitely seems to perform better
1522 // if (r) printf("R: %d\n",r);
1523 }
1524
1525 ds=H5Dcreate(file,ipath, hdt, spc, plist );
1526// printf("CRT %d %s\n",ds,ipath);
1527 H5Pclose(plist); // safe to do this here?
1528 }
1529 else { //existing file
1530 hid_t spc_file = H5Dget_space(ds);
1531 rank = H5Sget_simple_extent_ndims(spc_file);
1532 H5Sclose(spc_file);
1533 }
1534
1535 if (! data) {
1536 H5Dclose(ds);
1537 H5Sclose(spc);
1538// std::cerr << "Warning:blank image written!!! " << std::endl;
1539 return 0;
1540 }
1541
1542 // convert data to unsigned short, unsigned char...
1543 hsize_t size = (hsize_t)nx*ny*nz;
1544
1545 unsigned char *ucdata = NULL;
1546 unsigned short *usdata = NULL;
1547 char *cdata = NULL;
1548 short *sdata = NULL;
1549
1550 EMUtil::getRenderMinMax(data, nx, ny, rendermin, rendermax, renderbits, nz);
1551 int rendertrunc = 0; // keep track of truncated pixels
1552// printf("RMM %f %f\n",rendermin,rendermax);
1553
1554 // Limiting values for signed and unsigned with specified bits
1555 float RUMIN = 0;
1556 float RUMAX = (1<<renderbits)-1.0f;
1557 float RSMIN = -(1<<(renderbits-1));
1558 float RSMAX = (1<<(renderbits-1))-1;
1559
1560 bool scaled=0; // set if the data will need rescaling upon read
1561 herr_t err_no;
1562 if (area) {
1563 hid_t filespace = H5Dget_space(ds);
1564 hid_t memoryspace = 0;
1565
1566 if (rank == 3) {
1567 hsize_t doffset[3]; /* hyperslab offset in the file */
1568 doffset[0] = (hsize_t)(area->z_origin()<0 ? 0 : area->z_origin());
1569 doffset[1] = (hsize_t)(area->y_origin()<0 ? 0 : area->y_origin());
1570 doffset[2] = (hsize_t)(area->x_origin()<0 ? 0 : area->x_origin());
1571
1572 hsize_t dcount[3]; /* size of the hyperslab in the file */
1573 dcount[0] = (hsize_t)(area->get_depth()?area->get_depth():1);
1574 dcount[1] = (hsize_t)(area->get_height()?area->get_height():1);
1575 dcount[2] = (hsize_t)(area->get_width()?area->get_width():1);
1576
1577 herr_t err_no = H5Sselect_hyperslab(filespace, H5S_SELECT_SET,
1578 (const hsize_t*)doffset, NULL, (const hsize_t*)dcount, NULL);
1579
1580 if (err_no < 0) {
1581 std::cerr << "H5Sselect_hyperslab error: " << err_no << std::endl;
1582 }
1583
1584 /* Create memory space with size of the region. */
1585 hsize_t dims[3]; /* size of the region in the memory */
1586 dims[0] = dcount[2]?dcount[2]:1;
1587 dims[1] = dcount[1]?dcount[1]:1;
1588 dims[2] = dcount[0]?dcount[0]:1;
1589
1590 memoryspace = H5Screate_simple(rank, dims, NULL);
1591 }
1592 else if (rank == 2) {
1593 hsize_t doffset[2]; /* hyperslab offset in the file */
1594 doffset[0] = (hsize_t)(area->y_origin() < 0 ? 0 : area->y_origin());
1595 doffset[1] = (hsize_t)(area->x_origin() < 0 ? 0 : area->x_origin());
1596
1597 hsize_t dcount[2]; /* size of the hyperslab in the file */
1598 dcount[0] = (hsize_t)area->get_height();
1599 dcount[1] = (hsize_t)area->get_width();
1600
1601 herr_t err_no = H5Sselect_hyperslab(filespace, H5S_SELECT_SET,
1602 (const hsize_t*)doffset, NULL, (const hsize_t*)dcount, NULL);
1603
1604 if (err_no < 0) {
1605 std::cerr << "H5Sselect_hyperslab error: " << err_no << std::endl;
1606 }
1607
1608 /* Create memory space with size of the region. */
1609 /* Define memory dataspace - the memory we will created for the region */
1610 hsize_t dims[2]; /* size of the region in the memory */
1611 dims[0] = (hsize_t)(dcount[1]?dcount[1]:1);
1612 dims[1] = (hsize_t)(dcount[0]?dcount[0]:1);
1613
1614 memoryspace = H5Screate_simple(rank, dims, NULL);
1615 }
1616 else {
1617 std::cerr << "rank is wrong: " << rank << std::endl;
1618 }
1619
1620
1621 switch(dt) {
1622 case EMUtil::EM_FLOAT:
1623 err_no = H5Dwrite(ds, H5T_NATIVE_FLOAT, memoryspace, filespace, H5P_DEFAULT, data);
1624
1625 if (err_no < 0) {
1626 std::cerr << "H5Dwrite error float: " << err_no << std::endl;
1627 }
1628
1629 break;
1630 case EMUtil::EM_SHORT:
1631 sdata = new short[size];
1632
1633 for (size_t i = 0; i < size; ++i) {
1634 if (data[i] <= rendermin) {
1635 sdata[i] = (short)RSMIN;
1636 rendertrunc++;
1637 }
1638 else if (data[i] >= rendermax) {
1639 sdata[i] = (short)RSMAX;
1640 rendertrunc++;
1641 }
1642 else {
1643 sdata[i]=(short)roundf((data[i]-rendermin)/(rendermax-rendermin)*(RSMAX-RSMIN)+RSMIN);
1644 }
1645 }
1646
1647 err_no = H5Dwrite(ds, H5T_NATIVE_SHORT, memoryspace, filespace, H5P_DEFAULT, sdata);
1648
1649 if (err_no < 0) {
1650 std::cerr << "H5Dwrite error short: " << err_no << std::endl;
1651 }
1652
1653 if (sdata) {delete [] sdata; sdata = NULL;}
1654 scaled=1;
1655
1656 break;
1657 case EMUtil::EM_USHORT:
1658 usdata = new unsigned short[size];
1659
1660 for (size_t i = 0; i < size; ++i) {
1661 if (data[i] <= rendermin) {
1662 usdata[i] = 0;
1663 rendertrunc++;
1664 }
1665 else if (data[i] >= rendermax) {
1666 usdata[i] = (unsigned short)RUMAX;
1667 rendertrunc++;
1668 }
1669 else {
1670 usdata[i]=(unsigned short)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1671 }
1672 }
1673
1674 err_no = H5Dwrite(ds, H5T_NATIVE_USHORT, memoryspace, filespace, H5P_DEFAULT, usdata);
1675
1676 if (err_no < 0) {
1677 std::cerr << "H5Dwrite error ushort: " << err_no << std::endl;
1678 }
1679
1680 if (usdata) {delete [] usdata; usdata = NULL;}
1681 scaled=1;
1682
1683 break;
1684 case EMUtil::EM_CHAR:
1685 cdata = new char[size];
1686
1687 for (size_t i = 0; i < size; ++i) {
1688 if (data[i] <= rendermin) {
1689 cdata[i] = (char) RSMIN;
1690 rendertrunc++;
1691 }
1692 else if (data[i] >= rendermax){
1693 cdata[i] = (char) RSMAX;
1694 rendertrunc++;
1695 }
1696 else {
1697 cdata[i]=(char)roundf((data[i]-rendermin)/(rendermax-rendermin)*(RSMAX-RSMIN)+RSMIN);
1698 }
1699 }
1700
1701 err_no = H5Dwrite(ds, H5T_NATIVE_CHAR, memoryspace, filespace, H5P_DEFAULT, cdata);
1702
1703 if (err_no < 0) {
1704 std::cerr << "H5Dwrite error char: " << err_no << std::endl;
1705 }
1706
1707 if (cdata) {delete [] cdata; cdata = NULL;}
1708 scaled=1;
1709
1710 break;
1711 case EMUtil::EM_UCHAR:
1712 ucdata = new unsigned char[size];
1713
1714 for (size_t i = 0; i < size; ++i) {
1715 if (data[i] <= rendermin) {
1716 ucdata[i] = 0;
1717 rendertrunc++;
1718 }
1719 else if (data[i] >= rendermax){
1720 ucdata[i] = (unsigned char)RUMAX;
1721 rendertrunc++;
1722 }
1723 else {
1724 ucdata[i]=(unsigned char)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1725 }
1726 }
1727
1728 err_no = H5Dwrite(ds, H5T_NATIVE_UCHAR, memoryspace, filespace, H5P_DEFAULT, ucdata);
1729
1730 if (err_no < 0) {
1731 std::cerr << "H5Dwrite error uchar: " << err_no << std::endl;
1732 }
1733
1734 if (ucdata) {delete [] ucdata; ucdata = NULL;}
1735 scaled=1;
1736
1737 break;
1738 case EMUtil::EM_COMPRESSED:
1739 if (renderbits<=0) err_no = H5Dwrite(ds, H5T_NATIVE_FLOAT, memoryspace, filespace, H5P_DEFAULT, data);
1740 else if (renderbits<=8) {
1741 ucdata = new unsigned char[size];
1742
1743 for (size_t i = 0; i < size; ++i) {
1744 if (data[i] <= rendermin) {
1745 ucdata[i] = 0;
1746 rendertrunc++;
1747 }
1748 else if (data[i] >= rendermax){
1749 ucdata[i] = (unsigned char)RUMAX;
1750 rendertrunc++;
1751 }
1752 else {
1753 ucdata[i]=(unsigned char)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1754 }
1755 }
1756
1757 err_no = H5Dwrite(ds, H5T_NATIVE_UCHAR, memoryspace, filespace, H5P_DEFAULT, ucdata);
1758 scaled=1;
1759 if (ucdata) {delete [] ucdata; ucdata = NULL;}
1760 }
1761 else {
1762 usdata = new unsigned short[size];
1763
1764 for (size_t i = 0; i < size; ++i) {
1765 if (data[i] <= rendermin) {
1766 usdata[i] = 0;
1767 rendertrunc++;
1768 }
1769 else if (data[i] >= rendermax) {
1770 usdata[i] = (unsigned short)RUMAX;
1771 rendertrunc++;
1772 }
1773 else {
1774 usdata[i]=(unsigned short)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1775 }
1776 }
1777
1778 err_no = H5Dwrite(ds, H5T_NATIVE_USHORT, memoryspace, filespace, H5P_DEFAULT, usdata);
1779 scaled=1;
1780 if (usdata) {delete [] usdata; usdata = NULL;}
1781 }
1782
1783 if (err_no < 0) {
1784 std::cerr << "H5Dwrite error compressed: " << err_no << std::endl;
1785 }
1786 break;
1787 default:
1788 throw ImageWriteException(filename,"HDF5 does not support region writing for this data format");
1789 }
1790
1791 H5Sclose(filespace);
1792 H5Sclose(memoryspace);
1793 }
1794 else {
1795 switch(dt) {
1796 case EMUtil::EM_FLOAT:
1797 H5Dwrite(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
1798 break;
1799 case EMUtil::EM_SHORT:
1800 sdata = new short[size];
1801
1802 for (size_t i = 0; i < size; ++i) {
1803 if (data[i] <= rendermin) {
1804 sdata[i] = (short)RSMIN;
1805 rendertrunc++;
1806 }
1807 else if (data[i] >= rendermax) {
1808 sdata[i] = (short)RSMAX;
1809 rendertrunc++;
1810 }
1811 else {
1812 sdata[i]=(short)roundf((data[i]-rendermin)/(rendermax-rendermin)*(RSMAX-RSMIN)+RSMIN);
1813 }
1814 }
1815
1816 H5Dwrite(ds,H5T_NATIVE_SHORT,spc,spc,H5P_DEFAULT,sdata);
1817
1818 if (sdata) {delete [] sdata; sdata = NULL;}
1819 scaled=1;
1820
1821 break;
1822 case EMUtil::EM_USHORT:
1823 usdata = new unsigned short[size];
1824
1825 for (size_t i = 0; i < size; ++i) {
1826 if (data[i] <= rendermin) {
1827 usdata[i] = 0;
1828 rendertrunc++;
1829 }
1830 else if (data[i] >= rendermax) {
1831 usdata[i] = (unsigned short)RUMAX;
1832 rendertrunc++;
1833 }
1834 else {
1835 usdata[i]=(unsigned short)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1836 }
1837 }
1838
1839 H5Dwrite(ds,H5T_NATIVE_USHORT,spc,spc,H5P_DEFAULT,usdata);
1840
1841 if (usdata) {delete [] usdata; usdata = NULL;}
1842 scaled=1;
1843
1844 break;
1845 case EMUtil::EM_CHAR:
1846 cdata = new char[size];
1847
1848 for (size_t i = 0; i < size; ++i) {
1849 if (data[i] <= rendermin) {
1850 cdata[i] = (char) RSMIN;
1851 rendertrunc++;
1852 }
1853 else if (data[i] >= rendermax){
1854 cdata[i] = (char) RSMAX;
1855 rendertrunc++;
1856 }
1857 else {
1858 cdata[i]=(char)roundf((data[i]-rendermin)/(rendermax-rendermin)*(RSMAX-RSMIN)+RSMIN);
1859 }
1860 }
1861
1862 H5Dwrite(ds,H5T_NATIVE_CHAR,spc,spc,H5P_DEFAULT,cdata);
1863
1864 if (cdata) {delete [] cdata; cdata = NULL;}
1865 scaled=1;
1866
1867 break;
1868 case EMUtil::EM_UCHAR:
1869 ucdata = new unsigned char[size];
1870
1871 for (size_t i = 0; i < size; ++i) {
1872 if (data[i] <= rendermin) {
1873 ucdata[i] = 0;
1874 rendertrunc++;
1875 }
1876 else if (data[i] >= rendermax){
1877 ucdata[i] = (unsigned char)RUMAX;
1878 rendertrunc++;
1879 }
1880 else {
1881 ucdata[i]=(unsigned char)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1882 }
1883 }
1884
1885 H5Dwrite(ds,H5T_NATIVE_UCHAR,spc,spc,H5P_DEFAULT,ucdata);
1886
1887 if (ucdata) {delete [] ucdata; ucdata = NULL;}
1888 scaled=1;
1889
1890 break;
1891 case EMUtil::EM_COMPRESSED:
1892 //printf("writec %d %f %f\n",renderbits,rendermin,rendermax);
1893 if (renderbits<=0) err_no = H5Dwrite(ds,H5T_NATIVE_FLOAT,spc,spc,H5P_DEFAULT,data);
1894 else if (renderbits<=8) {
1895 ucdata = new unsigned char[size];
1896
1897 for (size_t i = 0; i < size; ++i) {
1898 if (data[i] <= rendermin) {
1899 ucdata[i] = 0;
1900 rendertrunc++;
1901 }
1902 else if (data[i] >= rendermax) {
1903 ucdata[i] = (unsigned char)RUMAX;
1904 rendertrunc++;
1905 }
1906 else {
1907 ucdata[i]=(unsigned char)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1908 }
1909 }
1910
1911 err_no = H5Dwrite(ds,H5T_NATIVE_UCHAR,spc,spc,H5P_DEFAULT,ucdata);
1912 scaled=1;
1913 if (ucdata) {delete [] ucdata; ucdata = NULL;}
1914 }
1915 else {
1916 usdata = new unsigned short[size];
1917
1918 for (size_t i = 0; i < size; ++i) {
1919 if (data[i] <= rendermin) {
1920 usdata[i] = 0;
1921 rendertrunc++;
1922 }
1923 else if (data[i] >= rendermax) {
1924 usdata[i] = (unsigned short)RUMAX;
1925 rendertrunc++;
1926 }
1927 else {
1928 usdata[i]=(unsigned short)roundf((data[i]-rendermin)/(rendermax-rendermin)*RUMAX);
1929 }
1930 }
1931
1932 err_no = H5Dwrite(ds,H5T_NATIVE_USHORT,spc,spc,H5P_DEFAULT,usdata);
1933 scaled=1;
1934 if (usdata) {delete [] usdata; usdata = NULL;}
1935 }
1936
1937 if (err_no < 0) {
1938 printf("%d %f %f\n",renderbits,rendermin,rendermax);
1939 std::cerr << "H5Dwrite error compressed full: " << err_no << std::endl;
1940 }
1941
1942 break;
1943 default:
1944 throw ImageWriteException(filename,"HDF5 does not support this data format");
1945 }
1946 }
1947
1948 H5Sclose(spc);
1949 H5Dclose(ds);
1950
1951 if (scaled) {
1952 sprintf(ipath,"/MDF/images/%d",image_index);
1953 hid_t igrp=H5Gopen(file,ipath);
1954 write_attr(igrp,"EMAN.stored_rendermin",EMObject(rendermin));
1955 write_attr(igrp,"EMAN.stored_rendermax",EMObject(rendermax));
1956 write_attr(igrp,"EMAN.stored_renderbits",EMObject(renderbits));
1957 write_attr(igrp,"EMAN.stored_truncated",EMObject(rendertrunc));
1958 H5Gclose(igrp);
1959 }
1960
1961
1962 EXITFUNC;
1963 return 0;
1964}
1965
1966int HdfIO2::get_nimg()
1967{
1968 init();
1969 hid_t attr=H5Aopen_name(group,"imageid_max");
1970 int n = read_attr(attr);
1971 H5Aclose(attr);
1972
1973 return n+1;
1974}
1975
1976void HdfIO2::flush()
1977{
1978 return;
1979}
1980
1981bool HdfIO2::is_complex_mode()
1982{
1983 return false;
1984}
1985
1986// always big endian
1987bool HdfIO2::is_image_big_endian()
1988{
1989 return true;
1990}
1991
1992#endif //USE_HDF5
EMAN::Ctf
Ctf is the base class for all CTF model.
Definition: ctf.h:60
EMAN::Dict
Dict is a dictionary to store <string, EMObject> pair.
Definition: emobject.h:385
EMAN::Dict::values
vector< EMObject > values() const
Get a vector containing copies of each of the EMObjects in this dictionary.
Definition: emobject.h:487
EMAN::Dict::has_key
bool has_key(const string &key) const
Ask the Dictionary if it as a particular key.
Definition: emobject.h:511
EMAN::Dict::keys
vector< string > keys() const
Get a vector containing all of the (string) keys in this dictionary.
Definition: emobject.h:475
EMAN::EMObject
EMObject is a wrapper class for types including int, float, double, etc as defined in ObjectType.
Definition: emobject.h:123
EMAN::EMObject::get_type
ObjectType get_type() const
Get the ObjectType This is an enumerated type first declared in the class EMObjectTypes.
Definition: emobject.h:224
EMAN::EMUtil::EMDataType
EMDataType
Image pixel data type used in EMAN.
Definition: emutil.h:92
EMAN::ImageIO
ImageIO classes are designed for reading/writing various electron micrography image formats,...
Definition: imageio.h:127
EMAN::IntSize
IntSize is used to describe a 1D, 2D or 3D rectangular size in integers.
Definition: geometry.h:49
EMAN::Region
Region defines a 2D or 3D rectangular region specified by its origin coordinates and all edges' sizes...
Definition: geometry.h:497
EMAN::Region::get_width
float get_width() const
get the width
Definition: geometry.h:606
EMAN::Region::y_origin
float y_origin() const
get the y element of the origin
Definition: geometry.h:622
EMAN::Region::x_origin
float x_origin() const
get the x element of the origin
Definition: geometry.h:620
EMAN::Region::get_depth
float get_depth() const
get the depth
Definition: geometry.h:610
EMAN::Region::origin
FloatPoint origin
Definition: geometry.h:654
EMAN::Region::z_origin
float z_origin() const
get the z element of the origin
Definition: geometry.h:624
EMAN::Region::get_height
float get_height() const
get the height
Definition: geometry.h:608
EMAN::Transform
A Transform object is a somewhat specialized object designed specifically for EMAN2/Sparx storage of ...
Definition: transform.h:75
EMAN::Transform::at
float at(int r, int c) const
Get the value stored in the internal transformation matrix at at coordinate (r,c)
Definition: transform.h:396
read_data
void read_data(string fsp, size_t loc, const Region *area=0, const int file_nx=0, const int file_ny=0, const int file_nz=0)
Read the image pixel data in native byte order from a disk file.
write_data
void write_data(string fsp, size_t loc, const Region *const area=0, const int file_nx=0, const int file_ny=0, const int file_nz=0)
Dump the image pixel data in native byte order to a disk file.
ATTR_NAME_LEN
static const int ATTR_NAME_LEN
Definition: emutil.cpp:56
ImageReadException
#define ImageReadException(filename, desc)
Definition: exception.h:204
FileAccessException
#define FileAccessException(filename)
Definition: exception.h:187
ImageWriteException
#define ImageWriteException(imagename, desc)
Definition: exception.h:223
LOGERR
#define LOGERR
Definition: log.h:51
ENTERFUNC
#define ENTERFUNC
Definition: log.h:48
EXITFUNC
#define EXITFUNC
Definition: log.h:49
EMAN
E2Exception class.
Definition: aligner.h:40
Generated by doxygen 1.9.2