ImageFactoryFromHDF5.h

 
#pragma once
 
#if defined(ImageFactoryFromHDF5_RECURSES)
#error Recursive header files inclusion detected in ImageFactoryFromHDF5.h
#else // defined(ImageFactoryFromHDF5_RECURSES)
#define ImageFactoryFromHDF5_RECURSES
 
#if !defined ImageFactoryFromHDF5_h
#define ImageFactoryFromHDF5_h
 
// Inclusions
#include <iostream>
#include "DGtal/base/Common.h"
#include "DGtal/base/ConceptUtils.h"
#include "DGtal/images/CImage.h"
#include "DGtal/base/Alias.h"
#include "DGtal/kernel/CBoundedNumber.h"
 
#include "hdf5.h"
 
namespace DGtal
{
 
  // template class H5DSpecializations
  template <typename TImageFactory, typename T>
  struct H5DSpecializations
  {
    // ----------------------- Standard services ------------------------------
 
    typedef TImageFactory ImageFactory;
    typedef typename ImageFactory::OutputImage::Value Value;
 
    static int H5DreadS(ImageFactory &anImageFactory, hid_t memspace, Value *data_out);
    static int H5DwriteS(ImageFactory &anImageFactory, hid_t memspace, Value *data_in);
 
  }; // end of class H5DSpecializations
 
  // template class H5DSpecializations
  template <typename TImageFactory>
  struct H5DSpecializations<TImageFactory, DGtal::uint8_t>
  {
    // ----------------------- Standard services ------------------------------
 
    typedef TImageFactory ImageFactory;
    typedef typename ImageFactory::OutputImage::Value Value;
 
    static int H5DreadS(ImageFactory &anImageFactory, hid_t memspace, Value *data_out)
    {
      return H5Dread(anImageFactory.dataset, H5T_NATIVE_UINT8, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_out);
    }
 
    static int H5DwriteS(ImageFactory &anImageFactory, hid_t memspace, Value *data_in)
    {
      return H5Dwrite(anImageFactory.dataset, H5T_NATIVE_UINT8, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_in);
    }
 
  }; // end of class H5DSpecializations
 
  // template class H5DSpecializations
  template <typename TImageFactory>
  struct H5DSpecializations<TImageFactory, DGtal::int32_t>
  {
    // ----------------------- Standard services ------------------------------
 
    typedef TImageFactory ImageFactory;
    typedef typename ImageFactory::OutputImage::Value Value;
 
    static int H5DreadS(ImageFactory &anImageFactory, hid_t memspace, Value *data_out)
    {
      return H5Dread(anImageFactory.dataset, H5T_NATIVE_INT32, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_out);
    }
 
    static int H5DwriteS(ImageFactory &anImageFactory, hid_t memspace, Value *data_in)
    {
      return H5Dwrite(anImageFactory.dataset, H5T_NATIVE_INT32, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_in);
    }
 
  }; // end of class H5DSpecializations
 
  // template class H5DSpecializations
  template <typename TImageFactory>
  struct H5DSpecializations<TImageFactory, DGtal::int64_t>
  {
    // ----------------------- Standard services ------------------------------
 
    typedef TImageFactory ImageFactory;
    typedef typename ImageFactory::OutputImage::Value Value;
 
    static int H5DreadS(ImageFactory &anImageFactory, hid_t memspace, Value *data_out)
    {
      return H5Dread(anImageFactory.dataset, H5T_NATIVE_INT64, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_out);
    }
 
    static int H5DwriteS(ImageFactory &anImageFactory, hid_t memspace, Value *data_in)
    {
      return H5Dwrite(anImageFactory.dataset, H5T_NATIVE_INT64, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_in);
    }
 
  }; // end of class H5DSpecializations
 
  // template class H5DSpecializations
  template <typename TImageFactory>
  struct H5DSpecializations<TImageFactory, double>
  {
    // ----------------------- Standard services ------------------------------
 
    typedef TImageFactory ImageFactory;
    typedef typename ImageFactory::OutputImage::Value Value;
 
    static int H5DreadS(ImageFactory &anImageFactory, hid_t memspace, Value *data_out)
    {
      return H5Dread(anImageFactory.dataset, H5T_NATIVE_DOUBLE, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_out);
    }
 
    static int H5DwriteS(ImageFactory &anImageFactory, hid_t memspace, Value *data_in)
    {
      return H5Dwrite(anImageFactory.dataset, H5T_NATIVE_DOUBLE, memspace, anImageFactory.dataspace, H5P_DEFAULT, data_in);
    }
 
  }; // end of class H5DSpecializations
 
  // Template class ImageFactoryFromHDF5
  template <typename TImageContainer>
  class ImageFactoryFromHDF5
  {
 
    // ----------------------- Types ------------------------------
 
  public:
    typedef ImageFactoryFromHDF5<TImageContainer> Self;
 
    BOOST_CONCEPT_ASSERT(( concepts::CImage<TImageContainer> ));
 
    typedef TImageContainer ImageContainer;
    typedef typename ImageContainer::Domain Domain;
 
    typedef ImageContainer OutputImage;
    typedef typename OutputImage::Value Value;
 
    BOOST_CONCEPT_ASSERT(( concepts::CBoundedNumber< Value > ));
 
    // ----------------------- Standard services ------------------------------
 
  public:
 
    ImageFactoryFromHDF5(const std::string & aFilename, const std::string & aDataset):
      myFilename(aFilename), myDataset(aDataset)
    {
      const int ddim = Domain::dimension;
 
      hsize_t dims_out[ddim];              // dataset dimensions
 
      // Open the file and the dataset.
      file = H5Fopen(myFilename.c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
      dataset = H5Dopen2(file, myDataset.c_str(), H5P_DEFAULT);
 
      // Get datatype and dataspace handles and then query dataset class, order, size, rank and dimensions.
      datatype = H5Dget_type(dataset); // datatype handle
 
      dataspace = H5Dget_space(dataset); // dataspace handle
 
      H5Sget_simple_extent_dims(dataspace, dims_out, NULL);
 
      // --
 
      typedef SpaceND<ddim> TSpace;
      typename TSpace::Point low, up;
 
      typename Domain::Integer d;
      for(d=0; d<ddim; d++)
      {
        low[d]=0;
        up[d]=dims_out[ddim-d-1]-1;
      }
 
      myDomain = new Domain(low, up);
    }
 
    ~ImageFactoryFromHDF5()
    {
      delete myDomain;
 
      // --
 
      // Close/release resources.
      H5Tclose(datatype);
      H5Dclose(dataset);
      H5Sclose(dataspace);
      H5Fclose(file);
    }
    
  private:
    
    ImageFactoryFromHDF5( const ImageFactoryFromHDF5 & other );
    
    ImageFactoryFromHDF5 & operator=( const ImageFactoryFromHDF5 & other );
 
    // ----------------------- Interface --------------------------------------
  public:
 
 
    const Domain & domain() const
    {
        return *myDomain;
    }
 
 
 
 
    void selfDisplay ( std::ostream & out ) const;
 
    bool isValid() const
    {
      return (myDomain->isValid());
    }
 
    OutputImage * requestImage(const Domain &aDomain) // time consuming
    {
      DGtal::IOException dgtalio;
      
      const int ddim = Domain::dimension;
 
      // --
 
      hsize_t offset[ddim];        // hyperslab offset in the file
      hsize_t count[ddim];         // size of the hyperslab in the file
 
      herr_t status;
      hsize_t dimsm[ddim];         // memory space dimensions
      hid_t memspace;
 
      hsize_t offset_out[ddim];    // hyperslab offset in memory
      hsize_t count_out[ddim];     // size of the hyperslab in memory
 
      int N_SUB[ddim];
      typename Domain::Integer d;
 
      int malloc_size=1;
      for(d=0; d<ddim; d++)
      {
        N_SUB[d] = (aDomain.upperBound()[ddim-d-1]-aDomain.lowerBound()[ddim-d-1])+1;
        malloc_size = malloc_size*N_SUB[d];
      }
 
      Value *data_out = (Value*) malloc (malloc_size * sizeof(Value)); // output buffer
      if (data_out == NULL)
      {
        trace.error() << "data_out malloc error in requestImage: " << (malloc_size * sizeof(Value)) << std::endl;
        throw dgtalio;
      }
 
      // Define hyperslab in the dataset.
      for(d=0; d<ddim; d++)
        offset[d] = aDomain.lowerBound()[ddim-d-1]-myDomain->lowerBound()[ddim-d-1];
      for(d=0; d<ddim; d++)
        count[d] = N_SUB[d];
      status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
      if (status)
      {
        trace.error() << " H5Sselect_hyperslab from dataspace error" << std::endl;
        throw dgtalio;
      }
 
      // Define the memory dataspace.
      for(d=0; d<ddim; d++)
        dimsm[d] = N_SUB[d];
      memspace = H5Screate_simple(ddim,dimsm,NULL);
 
      // Define memory hyperslab.
      for(d=0; d<ddim; d++)
        offset_out[d] = 0;
      for(d=0; d<ddim; d++)
        count_out[d] = N_SUB[d];
      status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, offset_out, NULL, count_out, NULL);
      if (status)
      {
        trace.error() << " H5Sselect_hyperslab from memspace error" << std::endl;
        throw dgtalio;
      }
 
      // Read data from hyperslab in the file into the hyperslab in memory.
      //status = H5Dread(dataset, H5T_NATIVE_INT, memspace, dataspace, H5P_DEFAULT, data_out);
      status = H5DSpecializations<Self, Value>::H5DreadS(*this, memspace, data_out);
      if (status)
      {
        trace.error() << " H5DSpecializations/H5DreadS error" << std::endl;
        throw dgtalio;
      }
 
      OutputImage* outputImage = new OutputImage(aDomain);
      if (outputImage == NULL)
      {
        trace.error() << "outputImage new error in requestImage: " << std::endl;
        throw dgtalio;
      }
 
      typedef SpaceND<ddim> TSpace;
      typename TSpace::Point a, b;
      for(d=0; d<ddim; d++)
      {
        a[d]=offset[ddim-d-1]+myDomain->lowerBound()[d];
        b[d]=a[d]+N_SUB[ddim-d-1]-1;
      }
      HyperRectDomain<TSpace> hrdomain(a,b);
 
      int p=0;
      for( typename HyperRectDomain<TSpace>::ConstIterator
            it = hrdomain.begin(), itend = hrdomain.end();
          it != itend;
          ++it)
      {
        outputImage->setValue((*it), data_out[ p++ ]);
      }
 
      H5Sclose(memspace);
 
      // --
 
      free(data_out);
 
      return outputImage;
    }
 
    void flushImage(OutputImage* outputImage)
    {
      DGtal::IOException dgtalio;
      
      const int ddim = Domain::dimension;
 
      // --
 
      hsize_t offset[ddim];        // hyperslab offset in the file
      hsize_t count[ddim];         // size of the hyperslab in the file
 
      herr_t status;
      hsize_t dimsm[ddim];         // memory space dimensions
      hid_t memspace;
 
      hsize_t offset_in[ddim];    // hyperslab offset in memory
      hsize_t count_in[ddim];     // size of the hyperslab in memory
 
      //int i[ddim];
      int N_SUB[ddim];
      typename Domain::Integer d;
 
      int malloc_size=1;
      for(d=0; d<ddim; d++)
      {
        N_SUB[d] = (outputImage->domain().upperBound()[ddim-d-1]-outputImage->domain().lowerBound()[ddim-d-1])+1;
        malloc_size = malloc_size*N_SUB[d];
      }
 
      Value *data_in = (Value*) malloc (malloc_size * sizeof(Value)); // input buffer
      if (data_in == NULL)
      {
        trace.error() << "data_in malloc error in flushImage: " << (malloc_size * sizeof(Value)) << std::endl;
        throw dgtalio;
      }
 
      // Define hyperslab in the dataset.
      for(d=0; d<ddim; d++)
        offset[d] = outputImage->domain().lowerBound()[ddim-d-1]-myDomain->lowerBound()[ddim-d-1];
      for(d=0; d<ddim; d++)
        count[d] = N_SUB[d];
      status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset, NULL, count, NULL);
      if (status)
      {
        trace.error() << " H5Sselect_hyperslab from dataspace error" << std::endl;
        throw dgtalio;
      }
 
      // Define the memory dataspace.
      for(d=0; d<ddim; d++)
        dimsm[d] = N_SUB[d];
      memspace = H5Screate_simple(ddim,dimsm,NULL);
 
      // Define memory hyperslab.
      for(d=0; d<ddim; d++)
        offset_in[d] = 0;
      for(d=0; d<ddim; d++)
        count_in[d] = N_SUB[d];
      status = H5Sselect_hyperslab(memspace, H5S_SELECT_SET, offset_in, NULL, count_in, NULL);
      if (status)
      {
        trace.error() << " H5Sselect_hyperslab from memspace error" << std::endl;
        throw dgtalio;
      }
 
      typedef SpaceND<ddim> TSpace;
      typename TSpace::Point a, b;
      for(d=0; d<ddim; d++)
      {
        a[d]=offset[ddim-d-1]+myDomain->lowerBound()[d];
        b[d]=a[d]+N_SUB[ddim-d-1]-1;
      }
      HyperRectDomain<TSpace> hrdomain(a,b);
 
      std::vector<typename TSpace::Dimension> v(ddim);
      for(d=0; d<ddim; d++)
        v[d]=d;
 
      int p=0;
      for( typename HyperRectDomain<TSpace>::ConstSubRange::ConstIterator
            it = hrdomain.subRange(v, a).begin(), itend = hrdomain.subRange(v, a).end();
          it != itend;
          ++it)
      {
        data_in[ p++ ] = outputImage->operator()((*it));
      }
 
      // Write data from hyperslab in memory into the hyperslab in the file.
      //status = H5Dwrite(dataset, H5T_NATIVE_INT, memspace, dataspace, H5P_DEFAULT, data_in);
      status = H5DSpecializations<Self, Value>::H5DwriteS(*this, memspace, data_in);
      if (status)
      {
        trace.error() << " H5DSpecializations/H5DwriteS error" << std::endl;
        throw dgtalio;
      }
 
      H5Sclose(memspace);
 
      // --
 
      free(data_in);
    }
 
    void detachImage(OutputImage* outputImage)
    {
      delete outputImage;
    }
 
    // ------------------------- Protected Datas ------------------------------
  private:
    //ImageFactoryFromHDF5() {}
 
    // ------------------------- Private Datas --------------------------------
  protected:
 
    Domain *myDomain;
 
    const std::string myFilename;
    const std::string myDataset;
 
  public:
 
    // HDF5 handles
    hid_t file, dataset;
    hid_t datatype, dataspace;
 
    // ------------------------- Internals ------------------------------------
  private:
 
  }; // end of class ImageFactoryFromHDF5
 
 
  template <typename TImageContainer>
  std::ostream&
  operator<< ( std::ostream & out, const ImageFactoryFromHDF5<TImageContainer> & object );
 
} // namespace DGtal
 
 
// Includes inline functions.
#include "DGtal/images/ImageFactoryFromHDF5.ih"
 
//                                                                           //
 
#endif // !defined ImageFactoryFromHDF5_h
 
#undef ImageFactoryFromHDF5_RECURSES
#endif // else defined(ImageFactoryFromHDF5_RECURSES)