Dans cet exemple, nous allons réaliser un programme qui produit des fichiers HDF5 (binaire structuré) contenant des tableaux à 2 dimensions de pattern spectra. Le programme prend 2 arguments :

• le nom du fichier contenant l'image en entrée qui construira un arbre MAX et des attributs
• le nom du fichier qui les données produites (principalement un pattern spectra)

Dans cet exemple nous proposons de réaliser un pattern spectra à 2 dimensions. Cela signifie que la contribution de chaque nœud sera cumulé en fonction

• de la taille du nœud
• et de sa caractéristique d'être compact ou non

Nous avons utilisé un script python pour visualiser le tableau à 2 dimensions. Nous avons joué avec les paramètres pour disposer d'axes linéaires ou logarithmiques.

	linear	log
linear
log

Il faut commencer par les déclarations d'usage. La première pour réaliser des traces lors de l'exécution (le calcul du périmètre et un peu lent). Les autres sont là pour rappel.

#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"

Nous aurons besoin des quelques classes fondamentales : la structure d'arbre, le poids des nœuds et un ensembles d'autres attributs pour le pattern spectra.

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

En fonction de ce qui a été déclaré, il peut être utile de simplifier les appels en omettant le contexte d'utilisation (par exemple en évitant de préfixer avec std::).

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

Dans notre exemple, nous nous limiterons à des images 2D et la construction d'un arbre MAX. Nous préfererons un histogramme logarithmique par défaut. Nous choisissons un pattern spectra de 900 casses (30 bin x 30 bin). La compression HDF5 a besoin de paramètres propres (taux de compression, taille de cache).

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

Nous déclarons la fonction patternSpectra qui réalisera le calcul et l'écriture du résultat. Elle sera donc dépendante de la nature des pixels lus InPixelT, de l'image en entrée et du nom de fichier à produire.

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

Nous définissons une fonction créant des intervalles (au nombre de binCount) de valeurs logarithmiques entre un minimum et un maximum.

/** 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;
}

Nous définissons une fonction créant des intervalles (au nombre de binCount) de valeurs de même taille entre un minimum et un maximum.

/** 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;
}

Le point d'entré du programme, commencera par analyser les arguments

• argv[0] : nom du programme
• argv[1] : nom du fichier contenant l'image en entrée
• argv[2] : nom du fichier contenant les données produites

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]);

En fonction des pixels de l'image nous invoquerons la fonction patternSpectra avec le bon type.

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

Voici la fonction principale du programme. Elle commence par lire les pixels.

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)
  const 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);

Nous déclarons l'objet arbre qui contiendra le chaînage des feuilles vers la racine.

  cout << "Build tree" << endl;
  // no border (i.e. all pixels are take in account)
  const Border<GeoDimT> border (size, false);
  // neighbors take in account (default connectivity C4)
  const 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));

Nous déclenchons la construction de l'arbre.

• ATB contient les éléments de l'images (matrice de pixels, connectivité, type d'arbre)
• tree coquille vide qui va contenir l'arbre
• weightAttributes est le premier tableau d'attributs associés aux nœuds (niveaux de gris)

  atb.buildTree (tree, weightAttributes);

Nous créons autant de tableau d'attributs associés aux nœuds que nécessaire (taille en pixels, périmètre, propriété d'être compact).

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

Nous construisons des intervalles de valeurs (absices et ordonnées).

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

Nous calculons le pattern spectra en parcourant tous les nœuds.

  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) {
    const DimBin areaBinId = lower_bound (areaScale.begin (), areaScale.end (), double (areaValues[compId])) - areaScale.begin ();
    const DimBin compactBinId = lower_bound (compactScale.begin (), compactScale.end (), double (compactValues[compId])) - compactScale.begin ();
    const double diff = abs (double (weightValues[compId]) -
			     double (weightValues[tree.getCompParent (compId)]));
    const double volume = diff * areaValues[compId];
 
    const size_t id (compactBinId + areaBinId * binCompactCount);
    psTab [id] += volume;
  }

Nous construisons, pour l'exemple, Une représentation réduite du pattern spectra (de nombreuses cases sont vides).

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

Nous sauvons les valeurs dans un fichier HDF5.

  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 contenues dans le fichier.

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

Le Pattern Spectra est constitué de 30 x 30 cases.

   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
      }
   }

Voici les information permettant de reconstituer l'axe des abscisses.

   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
      }
   }

Voici les information permettant de reconstituer l'axe des ordonnées.

   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
      }
   }

Pour l'exemple voici une autre représentation (réduite) du Parttern Spectra.

   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
      }
   }
}
}