DGtal  1.4.beta
HDF5Reader.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file HDF5Reader.ih
19  * @author Martial Tola (\c martial.tola@liris.cnrs.fr )
20  *
21  * @date 2013/04/16
22  *
23  * Implementation of inline methods defined in HDF5Reader.h
24  *
25  * This file is part of the DGtal library.
26  */
27 
28 ///////////////////////////////////////////////////////////////////////////////
29 // IMPLEMENTATION of inline methods.
30 ///////////////////////////////////////////////////////////////////////////////
31 
32 //////////////////////////////////////////////////////////////////////////////
33 #include <cstdlib>
34 #include <iostream>
35 #include <fstream>
36 #include <sstream>
37 
38 #include "DGtal/kernel/SpaceND.h"
39 #include "DGtal/kernel/domains/HyperRectDomain.h"
40 
41 #include <hdf5.h>
42 #include <hdf5_hl.h>
43 //////////////////////////////////////////////////////////////////////////////
44 
45 
46 
47 ///////////////////////////////////////////////////////////////////////////////
48 // Implementation of inline methods //
49 
50 
51 ///////////////////////////////////////////////////////////////////////////////
52 // Implementation of inline functions and external operators //
53 
54 template <typename TImageContainer, typename TFunctor>
55 inline
56 TImageContainer
57 DGtal::HDF5Reader<TImageContainer, TFunctor>::importHDF5(const std::string & aFilename, const std::string & aDataset,
58  const Functor & aFunctor, bool topbotomOrder)
59 {
60  DGtal::IOException dgtalio;
61 
62  BOOST_STATIC_ASSERT( (ImageContainer::Domain::dimension == 2));
63 
64  hid_t fileID;
65  herr_t hErrVal;
66  hsize_t width;
67  hsize_t height;
68  hsize_t planes;
69  char interlace[16];
70  hssize_t npals;
71  unsigned char *imgBuf; // 'unsigned char' is mandated by the HDF5 library
72 
73  // Load the library -- not required on most platforms
74  hErrVal = H5open();
75  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
76 
77  // Open an existing file
78  fileID = H5Fopen(aFilename.c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);
79  if (fileID < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
80 
81  // Figure out if the data is an image or not
82  hErrVal = H5IMis_image(fileID, aDataset.c_str());
83  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
84  if(hErrVal != 1) {
85  trace.error() << "The data is NOT an image" << std::endl;
86  throw dgtalio;
87  }
88 
89  // Get info about the image
90  interlace[0] = 0;
91  hErrVal = H5IMget_image_info (fileID, aDataset.c_str(), &width, &height, &planes, interlace, &npals);
92  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
93 
94  // If the image is the wrong type, then throw exception
95  if( !( ((planes==1) && (npals==1)) || ((planes==3) && (npals==0) && (strcmp(interlace, "INTERLACE_PIXEL")==0)) ) ) {
96  trace.error() << "ERROR: Image is not 8-bit with palette or image is not 24-bit truecolor with INTERLACE_PIXEL pixel layout" << std::endl;
97  throw dgtalio;
98  }
99 
100  // Allocate space for the image
101  imgBuf = (unsigned char*) malloc (width * height * planes * sizeof( unsigned char ));
102  if (imgBuf == NULL) {
103  trace.error() << "ERROR: Could not allocate space for image" << std::endl;
104  throw dgtalio;
105  }
106 
107  // Read the data from the file
108  hErrVal = H5IMread_image (fileID, aDataset.c_str(), imgBuf);
109  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
110 
111  // Close the file
112  hErrVal = H5Fclose(fileID);
113  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
114 
115  // Unload the library and free any remaining resources
116  hErrVal = H5close();
117  if (hErrVal < 0) { trace.error() << "HDF5 ERROR: " << std::endl; H5Eprint2(H5E_DEFAULT, stderr); throw dgtalio; }
118 
119  // --
120 
121  typename TImageContainer::Point firstPoint;
122  typename TImageContainer::Point lastPoint;
123 
124  firstPoint = TImageContainer::Point::zero;
125  lastPoint[0] = width-1;
126  lastPoint[1] = height-1;
127 
128  typename TImageContainer::Domain domain(firstPoint,lastPoint);
129  TImageContainer image(domain);
130 
131  for(hsize_t y=0; y<height; y++)
132  for(hsize_t x=0; x<width; x++)
133  {
134  typename TImageContainer::Point pt;
135  if (topbotomOrder){
136  pt[0]=x; pt[1]=height-1-y;
137  }else{
138  pt[0]=x; pt[1]=y;
139  }
140 
141  if( ((planes==1) && (npals==1)) )
142  image.setValue(pt, aFunctor(imgBuf[ y*width+x ]));
143  else if( ((planes==3) && (npals==0) && (strcmp(interlace, "INTERLACE_PIXEL")==0)) )
144  image.setValue(pt, aFunctor(imgBuf[ ( (y*planes*width+x*planes)+(y*planes*width+x*planes+1)+(y*planes*width+x*planes+2) )/3 ]));
145  }
146 
147  free(imgBuf);
148 
149  return image;
150 }
151 
152 template <typename TImageContainer, typename TFunctor>
153 inline
154 TImageContainer
155 DGtal::HDF5Reader<TImageContainer, TFunctor>::importHDF5_3D(const std::string & aFilename, const std::string & aDataset,
156  const Functor & aFunctor)
157 {
158  DGtal::IOException dgtalio;
159 
160  BOOST_STATIC_ASSERT( (ImageContainer::Domain::dimension == 3));
161 
162  const int ddim = ImageContainer::Domain::dimension;
163 
164  // image domain
165  Domain *myDomain;
166  Domain aDomain;
167 
168  // HDF5 handles
169  hid_t file, dataset;
170  hid_t datatype, dataspace;
171 
172  hsize_t dims_out[ddim]; // dataset dimensions
173 
174  // Open the file and the dataset.
175  file = H5Fopen(aFilename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
176  dataset = H5Dopen2(file, aDataset.c_str(), H5P_DEFAULT);
177 
178  // Get datatype and dataspace handles and then query dataset class, order, size, rank and dimensions.
179  datatype = H5Dget_type(dataset); // datatype handle
180 
181  dataspace = H5Dget_space(dataset); // dataspace handle
182 
183  H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
184 
185  // --
186 
187  typedef SpaceND<ddim> TSpace;
188  typename TSpace::Point low, up;
189 
190  typename Domain::Integer d;
191  for(d=0; d<ddim; d++)
192  {
193  low[d]=0;
194  up[d]=dims_out[ddim-d-1]-1;
195  }
196 
197  myDomain = new Domain(low, up);
198  aDomain = *myDomain; // because here we want the full image
199 
200  /* - reading HDF5- */
201  hsize_t offset[ddim]; // hyperslab offset in the file
202  hsize_t count[ddim]; // size of the hyperslab in the file
203 
204  herr_t status;
205  hsize_t dimsm[ddim]; // memory space dimensions
206  hid_t memspace;
207 
208  hsize_t offset_out[ddim]; // hyperslab offset in memory
209  hsize_t count_out[ddim]; // size of the hyperslab in memory
210 
211  //int i[ddim];
212  int N_SUB[ddim];
213 
214  int malloc_size=1;
215  for(d=0; d<ddim; d++)
216  {
217  N_SUB[d] = (aDomain.upperBound()[ddim-d-1]-aDomain.lowerBound()[ddim-d-1])+1;
218  malloc_size = malloc_size*N_SUB[d];
219  }
220 
221  Value *data_out = (Value*) malloc (malloc_size * sizeof(Value)); // output buffer
222  if (data_out == NULL)
223  {
224  trace.error() << "data_out malloc error in importHDF5_3D: " << (malloc_size * sizeof(Value)) << std::endl;
225  throw dgtalio;
226  }
227 
228  // Define hyperslab in the dataset.
229  for(d=0; d<ddim; d++)
230  offset[d] = aDomain.lowerBound()[ddim-d-1]-myDomain->lowerBound()[ddim-d-1];
231  for(d=0; d<ddim; d++)
232  count[d] = N_SUB[d];
233  status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
234  if (status)
235  {
236  trace.error() << " H5Sselect_hyperslab from dataspace error" << std::endl;
237  throw dgtalio;
238  }
239 
240  // Define the memory dataspace.
241  for(d=0; d<ddim; d++)
242  dimsm[d] = N_SUB[d];
243  memspace = H5Screate_simple(ddim,dimsm,NULL);
244 
245  // Define memory hyperslab.
246  for(d=0; d<ddim; d++)
247  offset_out[d] = 0;
248  for(d=0; d<ddim; d++)
249  count_out[d] = N_SUB[d];
250  status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, offset_out, NULL, count_out, NULL);
251  if (status)
252  {
253  trace.error() << " H5Sselect_hyperslab from memspace error" << std::endl;
254  throw dgtalio;
255  }
256 
257  // Read data from hyperslab in the file into the hyperslab in memory.
258  status = H5Dread(dataset, H5T_NATIVE_UINT8, memspace, dataspace, H5P_DEFAULT, data_out);
259  if (status)
260  {
261  trace.error() << " H5Dread error" << std::endl;
262  throw dgtalio;
263  }
264 
265  OutputImage outputImage(aDomain);
266 
267  typedef SpaceND<ddim> TSpace;
268  typename TSpace::Point a, b;
269  for(d=0; d<ddim; d++)
270  {
271  a[d]=offset[ddim-d-1]+myDomain->lowerBound()[d];
272  b[d]=a[d]+N_SUB[ddim-d-1]-1;
273  }
274  HyperRectDomain<TSpace> domain(a,b);
275 
276  int p=0;
277  for( typename HyperRectDomain<TSpace>/*::ConstSubRange*/::ConstIterator
278  it = domain/*.subRange(v, a)*/.begin(), itend = domain/*.subRange(v, a)*/.end();
279  it != itend;
280  ++it)
281  {
282  outputImage.setValue((*it), aFunctor(data_out[ p++ ]));
283  }
284 
285  H5Sclose(memspace);
286 
287  free(data_out);
288  /* - reading HDF5- */
289 
290  delete myDomain;
291 
292  // Close/release resources.
293  H5Tclose(datatype);
294  H5Dclose(dataset);
295  H5Sclose(dataspace);
296  H5Fclose(file);
297 
298  return outputImage;
299 }
300 
301 // //
302 ///////////////////////////////////////////////////////////////////////////////