Dans cet exemple, nous allons réaliser un programme qui produit des fichiers HDF5 (binaire structuré) contenant des tableau à 2 dimension de pattern spectra.

	linear	log
linear
log

#include <iostream>
 
// pixels type
#include "gdal.h"
// log functions
#include "obelixDebug.hpp"
// Point and Size (2D or 3D)
#include "obelixGeo.hpp"
// read/write images
#include "IImage.hpp"
// nodata pixels definition in images
#include "Border.hpp"
// tree structure
#include "Tree.hpp"

#include "ArrayTree/ArrayTreeBuilder.hpp"
#include "Attributes/WeightAttributes.hpp"
#include "Attributes/AreaAttributes.hpp"
#include "Attributes/PerimeterAttributes.hpp"
#include "Attributes/CompactnessAttributes.hpp"

using namespace std;
using namespace boost;
using namespace obelix;
using namespace obelix::triskele;

inline unsigned int abs (unsigned int x) { return x;}
 
typedef uint16_t DimBin;
 
static const int GeoDimT (2);
static TreeType treeType = MAX;
static DimBin linearCompactBin (false), linearAreaBin (false);
static DimBin binCompactCount (30), binAreaCount (30);
static unsigned int compressLevel (5);
static size_t cacheSize (1024*1024);

template<typename PixelT> void patternSpectra (const IImage<GeoDimT> &inputImage, const string &outputData);

/** Create log bin */
vector<double>
logSpace (const double &min, const double &max, const int &binCount) {
  vector<double> logSpace;
  logSpace.reserve (binCount);
  double exp_scale = (log (max) - log (min)) / (binCount-1);
  for (int i = 0; i < binCount; ++i)
    logSpace.push_back (exp (log (min) + i*exp_scale));
  return logSpace;
}

/** Create linear bin */
vector<double>
linearSpace (const double &min, const double &max, const int &binCount) {
  vector<double> linearSpace;
  linearSpace.reserve (binCount);
  for (int i = 0; i < binCount; ++i)
    linearSpace.push_back ( min + (max-min)*i/(binCount-1));
  return linearSpace;
}

int
main (int argc, char** argv, char** envp) {
  if (argc != 3) {
    cout << "Usage: " << argv[0] << " input-image-filename outpout-PS-filename" << endl;
    exit (1);
  }
  IImage<GeoDimT> inputImage (argv[1]);
  inputImage.readImage ();
  const string &outputData (argv[2]);

  switch (inputImage.getDataType ()) {
  case GDT_Byte:
    patternSpectra<uint8_t> (inputImage, outputData); break;
  case GDT_UInt16:
    patternSpectra<uint16_t> (inputImage, outputData); break;
  case GDT_Int16:
    patternSpectra<int16_t> (inputImage, outputData); break;
  case GDT_UInt32:
    patternSpectra<uint32_t> (inputImage, outputData); break;
  case GDT_Int32:
    patternSpectra<int32_t> (inputImage, outputData); break;
  case GDT_Float32:
    patternSpectra<float> (inputImage, outputData); break;
  case GDT_Float64:
    patternSpectra<double> (inputImage, outputData); break;
 
  default :
    cerr << "unknown type!" << endl; break;
    return 1;
  }
  return 0;
}

template<typename PixelT>
void
patternSpectra (const IImage<GeoDimT> &inputImage, const string &outputData) {
 
  cout << "Read image" << endl;
  // empty raster
  Raster<PixelT, GeoDimT> raster;
  // image size (width x height)
  Size<GeoDimT> size (inputImage.getSize ());
  // read first band (0) in 2D mode from origine [0,0] to end [width, height].
  inputImage.readBand (raster, 0, NullPoint2D, size);

  cout << "Build tree" << endl;
  // no border (i.e. all pixels are take in account)
  Border<GeoDimT> border (size, false);
  // neighbors take in account (default connectivity C4)
  GraphWalker<GeoDimT> graphWalker (border);
  // tree builder base on raster, connectivity and type of tree
  ArrayTreeBuilder<PixelT, PixelT, GeoDimT> atb (raster, graphWalker, treeType);
  // get the number of core
  const DimCore coreCount (boost::thread::hardware_concurrency ());
  // declare empty tree and set the number of thread for all algorithms
  Tree<GeoDimT> tree (coreCount);
  // reserve nodes weight attributes
  WeightAttributes<PixelT, GeoDimT> weightAttributes (tree, getDecrFromTreetype (treeType));

  atb.buildTree (tree, weightAttributes);

  // create area attribute
  cout << "Build Area Attributes" << endl;
  AreaAttributes<GeoDimT> areaAttributes (tree);
  cout << "Build Perimeter Attributes" << endl;
  PerimeterAttributes<GeoDimT> perimeterAttributes (tree, graphWalker);
  cout << "Build Compactness Attributes" << endl;
  CompactnessAttributes<GeoDimT> compactnessAttributes (tree, areaAttributes, perimeterAttributes);

  cout << "Build PS" << endl;
  static DimImg MaxArea (size.getPixelsCount ());
  vector<double> compactScale = linearCompactBin ?
    linearSpace (0., 1., binCompactCount) :
    logSpace (1./MaxArea, 1., binCompactCount);
  vector<double> areaScale = linearAreaBin ?
    linearSpace (1., MaxArea, binAreaCount) :
    logSpace (1, MaxArea, binAreaCount);

  vector<double> psTab (binAreaCount * binCompactCount, 0.);
  const DimParent rootId (tree.getRootId ());
  const DimImgPack *areaValues = areaAttributes.getValues ();
  const PixelT *weightValues = weightAttributes.getValues ();
  const double *compactValues = compactnessAttributes.getValues ();
 
  for (DimParent compId = 0; compId < rootId; ++compId) {
    DimBin areaBinId = lower_bound (areaScale.begin (), areaScale.end (), double (areaValues[compId])) - areaScale.begin ();
    DimBin compactBinId = lower_bound (compactScale.begin (), compactScale.end (), double (compactValues[compId])) - compactScale.begin ();
    double diff = abs (double (weightValues[compId])-
		       double (weightValues[tree.getCompParent (compId)]));
    double volume = diff * areaValues[compId];
 
    size_t id (compactBinId + areaBinId * binCompactCount);
    psTab [id] += volume;
  }

  vector<DimBin> reduceBinIdx;
  vector<double> reduceBinVal;
  for (DimBin binIdx = 0; binIdx < psTab.size (); ++binIdx) {
    if (!psTab[binIdx])
      continue;
    reduceBinIdx.push_back (binIdx);
    reduceBinVal.push_back (psTab[binIdx]);
  }

  cout << "Write HDF5" << endl;
  HDF5 hdf5;
  hdf5.set_compression_level (compressLevel);
  hdf5.set_cache_size (cacheSize);
  hdf5.open_write (outputData);
  hdf5.write_string ("Compactness", "XLabel", "/" );
  hdf5.write_string ("Area", "YLabel", "/" );
  hdf5.write_array (compactScale, "XTick", "/" );
  hdf5.write_array (areaScale, "Ytick", "/" );
  hdf5.write_array (psTab, "PS", "/" );
  hdf5.write_array (reduceBinIdx,"patchPSBinIdx","/" );
  hdf5.write_array (reduceBinVal,"patchPSBinVal","/" );
  hdf5.close_file ();
}

La commande h5ls permet de visualiser la structure du fichier

• PS : Tableau à 2 dimension contenant le Pattern Spectra (30*30 = 900)
• XLabel, XTick : Nom de l'axe des abscisses et des intervalles choisis
• YLabel, Ytick : Nom de l'axe des abscisses et des intervalles choisis
• patchPSBinIdx, patchPSBinVal : représentation réduite du tableau (seuls les cases avec une valeur sont enregistré avec leur index).

$ h5ls arles.kor 
PS                       Dataset {900}
XLabel                   Dataset {SCALAR}
XTick                    Dataset {30}
YLabel                   Dataset {SCALAR}
Ytick                    Dataset {30}
patchPSBinIdx            Dataset {69}
patchPSBinVal            Dataset {69}

La commande h5dump permet d'afficher les valeurs dans le fichier

$ h5dump arles.kor 
HDF5 "arles.kor" {
GROUP "/" {

Pattern Spectra

   DATASET "PS" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 900 ) / ( 900 ) }
      DATA {
      (0): 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
... 
      (883): 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
      }
   }

Abscisse

   DATASET "XLabel" {
      DATATYPE  H5T_STRING {
         STRSIZE 12;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "Compactness"
      }
   }
   DATASET "XTick" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 30 ) / ( 30 ) }
      DATA {
      (0): 0, 0.0344828, 0.0689655, 0.103448, 0.137931, 0.172414, 0.206897,
... 
      (25): 0.862069, 0.896552, 0.931034, 0.965517, 1
      }
   }

Ordonnée

   DATASET "YLabel" {
      DATATYPE  H5T_STRING {
         STRSIZE 5;
         STRPAD H5T_STR_NULLTERM;
         CSET H5T_CSET_ASCII;
         CTYPE H5T_C_S1;
      }
      DATASPACE  SCALAR
      DATA {
      (0): "Area"
      }
   }
   DATASET "Ytick" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 30 ) / ( 30 ) }
      DATA {
      (0): 1, 1.75326, 3.07394, 5.38942, 9.44908, 16.5667, 29.0459, 50.9251,
... 
      (28): 6.72556e+06, 1.17917e+07
      }
   }

Données réduites

   DATASET "patchPSBinIdx" {
      DATATYPE  H5T_STD_U16LE
      DATASPACE  SIMPLE { ( 69 ) / ( 69 ) }
      DATA {
      (0): 70, 75, 96, 98, 100, 124, 125, 126, 127, 152, 153, 154, 155, 156,
... 
      (66): 841, 871, 872
      }
   }
   DATASET "patchPSBinVal" {
      DATATYPE  H5T_IEEE_F64LE
      DATASPACE  SIMPLE { ( 69 ) / ( 69 ) }
      DATA {
      (0): 6.7979e+06, 7.63336e+06, 6.20837e+06, 7.33739e+06, 125310,
... 
      (66): 3.50345e+09, 2.9508e+09, 1.09653e+09
      }
   }
}
}