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 ();
}