====== Exemple Pattern Spectra ======
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
===== Visualisation =====
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 ^ {{ :cpp:arlestt.png?nolink&200 |}} | {{ :cpp:arlestf.png?nolink&200 |}} |
^ log ^ {{ :cpp:arlesft.png?nolink&200 |}} | {{ :cpp:arlesff.png?nolink&200 |}} |
===== Code =====
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
// 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 void patternSpectra (const IImage &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
logSpace (const double &min, const double &max, const int &binCount) {
vector 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
linearSpace (const double &min, const double &max, const int &binCount) {
vector 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ée 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 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 (inputImage, outputData); break;
case GDT_UInt16:
patternSpectra (inputImage, outputData); break;
case GDT_Int16:
patternSpectra (inputImage, outputData); break;
case GDT_UInt32:
patternSpectra (inputImage, outputData); break;
case GDT_Int32:
patternSpectra (inputImage, outputData); break;
case GDT_Float32:
patternSpectra (inputImage, outputData); break;
case GDT_Float64:
patternSpectra (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
void
patternSpectra (const IImage &inputImage, const string &outputData) {
cout << "Read image" << endl;
// empty raster
Raster raster;
// image size (width x height)
const Size 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 border (size, false);
// neighbors take in account (default connectivity C4)
const GraphWalker graphWalker (border);
// tree builder base on raster, connectivity and type of tree
ArrayTreeBuilder 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 tree (coreCount);
// reserve nodes weight attributes
WeightAttributes 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 areaAttributes (tree);
cout << "Build Perimeter Attributes" << endl;
const PerimeterAttributes perimeterAttributes (tree, graphWalker);
cout << "Build Compactness Attributes" << endl;
const CompactnessAttributes compactnessAttributes (tree, areaAttributes, perimeterAttributes);
Nous construisons des intervalles de valeurs (abscisses et ordonnées).
cout << "Build PS" << endl;
const DimImg MaxArea (size.getPixelsCount ());
const vector compactScale = linearCompactBin ?
linearSpace (0., 1., binCompactCount) :
logSpace (1./MaxArea, 1., binCompactCount);
const vector areaScale = linearAreaBin ?
linearSpace (1., MaxArea, binAreaCount) :
logSpace (1, MaxArea, binAreaCount);
Nous calculons le pattern spectra en parcourant tous les nœuds.
vector 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 reduceBinIdx;
vector 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 ();
}
===== Résultat =====
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
}
}
}
}