greyLevelBlob.h

Go to the documentation of this file.

/* This software and supporting documentation are distributed by
 *     Institut Federatif de Recherche 49
 *     CEA/NeuroSpin, Batiment 145,
 *     91191 Gif-sur-Yvette cedex
 *     France
 *
 * This software is governed by the CeCILL-B license under
 * French law and abiding by the rules of distribution of free software.
 * You can  use, modify and/or redistribute the software under the
 * terms of the CeCILL-B license as circulated by CEA, CNRS
 * and INRIA at the following URL "http://www.cecill.info".
 *
 * As a counterpart to the access to the source code and  rights to copy,
 * modify and redistribute granted by the license, users are provided only
 * with a limited warranty  and the software's author,  the holder of the
 * economic rights,  and the successive licensors  have only  limited
 * liability.
 *
 * In this respect, the user's attention is drawn to the risks associated
 * with loading,  using,  modifying and/or developing or reproducing the
 * software by the user in light of its specific status of free software,
 * that may mean  that it is complicated to manipulate,  and  that  also
 * therefore means  that it is reserved for developers  and  experienced
 * professionals having in-depth computer knowledge. Users are therefore
 * encouraged to load and test the software's suitability as regards their
 * requirements in conditions enabling the security of their systems and/or
 * data to be ensured and,  more generally, to use and operate it in the
 * same conditions as regards security.
 *
 * The fact that you are presently reading this means that you have had
 * knowledge of the CeCILL-B license and that you accept its terms.
 */
 
#ifndef AIMS_PRIMALSKETCH_GREYLEVELBLOB_H
#define AIMS_PRIMALSKETCH_GREYLEVELBLOB_H
 
#include <cstdlib>
#include <aims/config/aimsalgo_config.h>
#include <iostream>
#include <set>
#include <map>
#include <list>
#include <vector>
 
namespace aims {
    const int BACKGROUND = 0;
 
    template<typename T> struct ltstr_p3d
    { // comparaison de Sites pour conteneurs tri�s (set, map)
        bool operator()(const T p1, const T p2) const { return false;}
    };
 
    template <> struct ltstr_p3d<Point3d>
    {
        bool operator ()( const Point3d p1, const Point3d p2 ) const
        {
            if ( p1[2] < p2[2] )
                return true;
            else if ( p1[2] > p2[2] )
                return false;
            else {
                if ( p1[1] < p2[1] )
                    return true;
                else if ( p1[1] > p2[1] )
                    return false;
                else
                {
                    if ( p1[0] < p2[0] )
                        return true;
                    else return false;
                }
            }
        }
    };
 
    template <> struct ltstr_p3d<std::pair<Point3df, uint> >
    {
        bool operator () ( const std::pair<Point3df, uint> p1,
                           const std::pair<Point3df, uint> p2 ) const
        {
            if ( p1.second < p2.second )
                return true;
            else
                return false;
        }
    };
 
    template<class T> class MaximumPoint;
    template<class T> class SaddlePoint;
    template<class T> class GreyLevelBlob;
 
    class BlobMeasurements
    {
 
        public:
            BlobMeasurements() {}
            BlobMeasurements(const BlobMeasurements &other) {maxIntensity=other.maxIntensity; meanIntensity=other.meanIntensity; maxContrast=other.maxContrast; meanContrast= other.meanContrast;  area=other.area; t=other.t; tValue=other.tValue;}
            BlobMeasurements(float maxInt, float meanInt, float maxCont, float meanCont, float a, float tv=0.0, float tvalue=0.0):
                maxIntensity(maxInt), meanIntensity(meanInt), maxContrast(maxCont), meanContrast(meanCont), area(a), t(tv), tValue(tvalue) {}
 
            float maxIntensity;
            float meanIntensity;
            float maxContrast;
            float meanContrast;
            float area;
            float t;
            float tValue;
 
            BlobMeasurements & operator = (const BlobMeasurements & other);
    };

 
    template<class T> class GreyLevelBlob
    {
 
        protected:
 
            int             _label;
            std::set<T,ltstr_p3d<T> >          listePoints;
            SaddlePoint<T>  *saddle;
            MaximumPoint<T> *maximum;
            float           scale;
            bool            grow;
 
        public:
 
            GreyLevelBlob ( MaximumPoint<T> *node, int label )
            {
                saddle = 0;
                maximum = node;
                node->blob = this;
                grow = true;
                _label = label;
                AddPoint (maximum->_node);
            }
            GreyLevelBlob () {
                saddle = 0;
                maximum = 0;
                grow = true;
                _label = BACKGROUND;
            }
 
            BlobMeasurements    measurements;
 
            void AddPoint( T node ) { listePoints.insert(node); }
            void StopGrowth() { grow = false; }
            bool CanGrow() { return grow; }
            int  Label() { return _label; }
 
            std::set<T, ltstr_p3d<T> > & GetListePoints () { return listePoints; }
 
            void SetSaddle (SaddlePoint<T>  *node)
            {
                saddle = node;
                node->AddBlob (this);
                StopGrowth();
            }
 
            SaddlePoint<T> *GetSaddle () { return saddle; }
            MaximumPoint<T> *GetMaximum () { return maximum; }
            float GetScale () { return scale; }
            void SetScale (float t) { scale = t; }
 
            GreyLevelBlob<T> & operator = ( const GreyLevelBlob<T> & other );
 
    };

 
    template<class T> class SaddlePoint
    {
 
        public:
            T _node;
            std::list<GreyLevelBlob<T>* > blobList;
 
            SaddlePoint () { }
            SaddlePoint (T node) : _node(node) { }
            ~SaddlePoint () { }
 
            void AddBlob(GreyLevelBlob<T> *blob) { blobList.push_back(blob); }
            std::list<GreyLevelBlob<T>*> & GetBlobList() { return blobList; }
 
            SaddlePoint<T> & operator = ( const SaddlePoint<T>& other );
 
    };

 
    template<class T> class MaximumPoint
    {
 
        public:
            T                _node;
            GreyLevelBlob<T> *blob;
            ~MaximumPoint() { }
 
            MaximumPoint () : blob(0) { }
            MaximumPoint ( T node ) : _node(node), blob(0) { }
            MaximumPoint<T> & operator = ( const MaximumPoint<T>& other );
 
    };

 
  template<typename Geom, typename Text> class ExtractGreyLevelBlobs
  {
 
    private :
 
      typedef typename SiteType<Geom>::type Site;
      typedef typename TexType<Text>::type Val;
 
      TexturedData<Geom,Text> *_texdata;
      TexturedData<Geom,Text> *_rawtexdata;
      TexturedData<Geom, Text> *_mask;      // For convenience, data and mask are same type.
                                                       // This is the only simple solution to keep type genericity
                                                 // Mask is 0 outside, >0 inside
      std::map<int, GreyLevelBlob<Site>* >  blobMap; // map label -> blob
      std::list<SaddlePoint<Site> *>    saddleList;
      std::list<MaximumPoint<Site> *>   maximumList;
      TexturedData<Geom, Text> labelsImage; // mapping site->label (blob, background...)
 
      char *_stats;        // File name for blob statistics
 
      void DoMasking(); // if required, use a mask to remove some blobs
      void CheckMask(); // basic check : is the geometry of the mask the same than the data ?
 
    public :
 
      ExtractGreyLevelBlobs(TexturedData<Geom,Text> *texdata, TexturedData<Geom,Text> *rawtexdata, TexturedData<Geom,Text> *mask)
      : _texdata(texdata), _rawtexdata(rawtexdata), _mask(mask), _stats(0){ labelsImage=TexturedData<Geom, Text>(*_texdata); CheckMask();}
      ExtractGreyLevelBlobs(TexturedData<Geom,Text> *texdata, TexturedData<Geom,Text> *rawtexdata, TexturedData<Geom,Text> *mask, char *stats)
      : _texdata(texdata), _rawtexdata(rawtexdata), _mask(mask), _stats(stats) { labelsImage=TexturedData<Geom, Text>(*_texdata); CheckMask();}
 
 
      ExtractGreyLevelBlobs(TexturedData<Geom,Text> *texdata)
      : _texdata(texdata), _rawtexdata(NULL), _mask(0), _stats(0) {                    labelsImage=TexturedData<Geom, Text>(*_texdata);}
      ExtractGreyLevelBlobs(TexturedData<Geom,Text> *texdata, TexturedData<Geom,Text> *mask)
      : _texdata(texdata), _rawtexdata(NULL),  _mask(mask), _stats(0) { labelsImage=TexturedData<Geom, Text>(*_texdata); CheckMask();}
      ExtractGreyLevelBlobs(TexturedData<Geom,Text> *texdata, TexturedData<Geom,Text> *mask, char *stats)
      : _texdata(texdata), _rawtexdata(NULL),  _mask(mask), _stats(stats) { labelsImage=TexturedData<Geom, Text>(*_texdata); CheckMask();}
      void Run(); // possibility of masking to exclude an area
 
       void ComputeBlobMeasurements();
       void SetOriginalTexture ( TexturedData<Geom,Text> *rawtexdata ) { _rawtexdata = rawtexdata; }
      TexturedData<Geom,Text> *GetOriginalTexture ( ) { return _rawtexdata; }
    // doit renvoyer une liste de GLB, de maximums et de saddle points
      std::list<SaddlePoint<Site> *> GetSaddleList() { return saddleList; }
      std::list<MaximumPoint<Site> *> GetMaxList() { return  maximumList; }
      std::map<int, GreyLevelBlob<Site> *> GetBlobs() { return blobMap; }
      TexturedData<Geom, Text> GetBlobImage() { return labelsImage; }
      int nbBlobs() { return blobMap.size(); }
 
  };
 

 
 template<typename T> class GreyLevelBlobTools
 {
 
     private :
     GreyLevelBlob<T> *_blob;
 
     public :
 
     GreyLevelBlobTools(GreyLevelBlob<T> *blob) : _blob(blob) {}
     std::vector<float> Boundingbox();
     std::vector<float> Barycenter();
 
 };
 
  //---------------------------------------------------------------------------
  //
  // Methods
  //
  //---------------------------------------------------------------------------
 
    inline
    BlobMeasurements& BlobMeasurements::operator = (const BlobMeasurements & other)
    {
        maxIntensity=other.maxIntensity;
        meanIntensity=other.meanIntensity;
        maxContrast=other.maxContrast;
        meanContrast=other.meanContrast;
        area=other.area;
        t=other.t;
        tValue=other.tValue;
 
     return *this;
    }
 
      //----------------------------------------------------------------------------------------------------------
 
     template<class T>
     inline
     SaddlePoint<T>& SaddlePoint<T>::operator = (const SaddlePoint<T> & other)
     {
          _node=other.node;
          return *this;
     }
 
       //----------------------------------------------------------------------------------------------------------
 
     template<class T>
     inline
     MaximumPoint<T>& MaximumPoint<T>::operator = (const MaximumPoint<T> & other)
     {
          _node=other.node;
          blob=other.blob;
          return *this;
     }
 
       //----------------------------------------------------------------------------------------------------------
 
     template<class T>
     inline
     GreyLevelBlob<T>& GreyLevelBlob<T>::operator = (const GreyLevelBlob<T> & other)
     {
          _label=other._label;
          listePoints=other.listePoints;
          scale=other.scale;
          saddle=other.saddle;
          maximum=other.maximum;
 
          return *this;
     }
 
       //----------------------------------------------------------------------------------------------------------
 
    template<typename Geom, typename Text> void ExtractGreyLevelBlobs<Geom, Text>::Run()
    {
        int currentLabel=0, labelBlob;
        MaximumPoint<Site> *maximum;
        SaddlePoint<Site> *saddle;
        GreyLevelBlob<Site> *greyLevelBlob;
 
        // multimap<intensity, data.site> for intensity based sorting of sites/voxels
        std::multimap<const Val, Site> sortedSites;
        SiteIterator<Geom> ptSite=(*_texdata).siteBegin();
 
        for ( ; ptSite != (*_texdata).siteEnd(); ++ptSite)
            sortedSites.insert(std::pair<const Val, Site>((*_texdata).intensity(*ptSite), (*ptSite)));
 
 
 
        // TexturedData that will contains the label associated with each site,
        // that is:  BACKGROUND, or grey-level blob number.
        // This TexturedData also gives acess to the neighbours of each site
 
        // Rem: labels should be (int) but will be same type than Text to handle both surfaces
        // and images. Not elegant but that's all I found so far.
 
        for (ptSite=labelsImage.siteBegin(); ptSite !=  labelsImage.siteEnd(); ++ptSite)
            labelsImage.intensity(*ptSite)=BACKGROUND;
        // GL Blob extraction
        // This is the core algorithm
        typename std::multimap<const Val, Site>::reverse_iterator ptSortedSites;
 
        //std::cout << "DEBUG: extractor has started detection" << std::endl;
        for ( ptSortedSites=sortedSites.rbegin(); ptSortedSites != sortedSites.rend(); ++ptSortedSites) {
            Val intensity=(*_texdata).intensity((*ptSortedSites).second);
 
            int nbAbove;
            int nbAboveLabels;
 
            // How many neighbours with higher intensity ?
 
            std::vector<Site> neighb=(*_texdata).neighbours((*ptSortedSites).second); //neighbours
 
            std::vector<Site> above; //neighbours with higher intensity
            std::set<int> aboveLabels; // labels of (above)
            typename std::vector<Site>::iterator ptNeigh=neighb.begin();
 
            for ( ; ptNeigh != neighb.end(); ++ptNeigh){
 
                if ((*_texdata).intensity(*ptNeigh) > intensity){
 
                    above.push_back(*ptNeigh);
                    aboveLabels.insert(int(labelsImage.intensity(*ptNeigh)));
                }
            }
 
            nbAbove=above.size();
            nbAboveLabels=aboveLabels.size();
 
            if (nbAbove==0){      // Case 1: local max, new blob
                currentLabel++;
                maximum=new MaximumPoint<Site>((*ptSortedSites).second);
                greyLevelBlob=new GreyLevelBlob<Site>(maximum, currentLabel);
                labelsImage.intensity((*ptSortedSites).second)=(Val) currentLabel;
                maximumList.push_back(maximum);
                blobMap[currentLabel]=(greyLevelBlob);
 
            }
            else if ( (nbAbove==1) && ((labelsImage.intensity(*(above.begin())))==BACKGROUND) ){
                                // Case 2: point in background
                labelsImage.intensity((*ptSortedSites).second)=BACKGROUND;
 
            }
 
            else if ((nbAbove>1) && (nbAboveLabels>1)){
                               // Case 3: one or several blobs to stop -> Saddle point
 
                int flagS=0;
                labelsImage.intensity((*ptSortedSites).second)=BACKGROUND;
                std::set<int>::iterator ptLabels=aboveLabels.begin();
                for (; ptLabels != aboveLabels.end(); ++ptLabels)
                if (*ptLabels!=BACKGROUND)
                    if (blobMap[*ptLabels]->CanGrow())
                        {
                            if (flagS==0)
                            {
                        saddle=new SaddlePoint<Site>((*ptSortedSites).second);
                        saddleList.push_back(saddle);
                            blobMap[*ptLabels]->SetSaddle(saddle);
                                flagS=1;
                            }
                            else
                                blobMap[*ptLabels]->SetSaddle(saddle);
                        }
            }
            else                 // Case 4: point belong to a growing blob
            {
                labelBlob=int(labelsImage.intensity(*(above.begin())));
                if ((labelBlob != BACKGROUND) && (blobMap[labelBlob]->CanGrow() ))
                {
                    labelsImage.intensity((*ptSortedSites).second)=(Val) labelBlob;
                    blobMap[labelBlob]->AddPoint((*ptSortedSites).second);
                }
                else
                {
                    labelsImage.intensity((*ptSortedSites).second)=BACKGROUND;
                }
            }
        }
        //std::cout << "DEBUG: found " << blobMap.size() << " blobs" << std::endl;
        //std::cout << "DEBUG: extractor about to enter masking" << std::endl;
        if (_mask!=0) { DoMasking();
            std::cout<<"MASKING" << std::endl;
        }
        //std::cout << "DEBUG: masking done, extractor about to enter measurements" << std::endl;
        ComputeBlobMeasurements();
        //std::cout << "DEBUG: done" << std::endl;
    }
 
    //----------------------------------------------------------------------------------------------------------
 
    template<typename Geom, typename Text> void ExtractGreyLevelBlobs<Geom, Text>::DoMasking() {
        typename std::list<MaximumPoint<Site> *>::iterator itMax=maximumList.begin(), itPrevious=itMax;
        SaddlePoint<Site> *saddle;
        double EPSILON = 0.00001;
        int del = 0;
    
        std::cout << "DEBUG: found " << blobMap.size() << " blobs" << std::endl;
        std::cout << "DEBUG: saddle list size : " << saddleList.size() << std::endl;
    
        for ( ; itMax != maximumList.end() ; ++itMax ) {
            if ( del == 1 ) {
                maximumList.erase(itPrevious);
            }
            del = 0;
            if ( fabs ( _mask->intensity((*itMax)->_node)) < EPSILON ) { // this is to test float=0.0             
                blobMap.erase ( (*itMax)->blob->Label() );                // remove the blob from the list.
    
                saddle = (*itMax)->blob->GetSaddle();                    // check the saddle
                if ( saddle != 0 ) { // some blobs can have no saddle if they are on the image border
                    if ( saddle->blobList.size() > 1 ) {
                        saddle->blobList.remove( (*itMax)->blob );
                    }
                    else {
                        saddle->blobList.remove((*itMax)->blob);
                        saddleList.remove(saddle);
                        delete saddle;
                    } 
                }
                delete (*itMax)->blob;        // delete the blob
                del = 1;                             // maximum has to be deleted
            }
            itPrevious=itMax;
        }
        if ( del == 1 )
            maximumList.erase(itPrevious);
        std::cout << "DEBUG: found " << blobMap.size() << " blobs" << std::endl;
    }
 
    //----------------------------------------------------------------------------------------------------------
 
  template<typename Geom, typename Text> void ExtractGreyLevelBlobs<Geom, Text>::CheckMask()
  {
          if (_mask != 0)
               if (_mask->NbSites() != _texdata->NbSites())
               {
                    std::cerr << "ExtractGreyLevelBlobs : mask and data do not have the same number of sites !!!" << std::endl;
                    exit(EXIT_FAILURE);
               }
  }
 
  //----------------------------------------------------------------------------------------------------------
 
  template<typename Geom, typename Text>
  void
  ExtractGreyLevelBlobs<Geom, Text>::ComputeBlobMeasurements ( ) {
      typename std::map<int, GreyLevelBlob<Site>* >::iterator blobIt;
      GreyLevelBlob<Site> *blob;
      SaddlePoint<Site> *saddle;
      MaximumPoint<Site> *maxi;
      Val saddleInt, maxInt, intensity;
      float maxIntensity, meanIntensity, maxContrast, meanContrast, area;
      int nbPoint;
      FILE *fileStat = 0;
 
      if ( _stats != 0 ) {
          fileStat = fopen ( _stats, "w" );
          fprintf ( fileStat, "Blob_label max_int mean_int max_cont mean_cont area\n" );
      }
 
      //std::cout << "DEBUG: going through blob list" << std::endl;
      //std::cout << "DEBUG: Size of list : " << blobMap.size() << std::endl;
      //std::cout << "DEBUG: saddle list size : " << saddleList.size() << std::endl;
      for ( blobIt = blobMap.begin() ; blobIt != blobMap.end() ; ++blobIt ) {
          //std::cout << ((*blobIt).second)->Label(); fflulongsh(stdout);
          //std::cout << "DEBUG : blob " << (*blobIt).first << std::endl;
          maxIntensity = 0.0;
          meanIntensity = 0.0;
          maxContrast = 0.0;
          meanContrast = 0.0;
          nbPoint = 0;
          area = 0.0; // pour l'instant il s'agit du nombre de points
                              // donc tres mauvais pour les surfaces
          blob = (*blobIt).second;
          //std::cout << "DEBUG : getting saddle" << std::endl;
          saddle = blob->GetSaddle();
          //std::cout << "DEBUG : getting maxi" << std::endl;
          maxi=blob->GetMaximum();
          //std::cout << "DEBUG : sites : " << (*_texdata).NbSites() << std::endl;
          //std::cout << "DEBUG : getting max intensity" << std::endl;
          maxInt=(*_texdata).intensity(maxi->_node);
          //std::cout << "DEBUG : getting saddle intensity" << std::endl;
          if (saddle != 0) {
              //std::cout << "DEBUG : saddle : " << (saddle) << std::endl;
              saddleInt=(*_texdata).intensity(saddle->_node);
          }
          else {
                //std::cout << "DEBUG : blob " << blob->Label() << " has no saddle " << std::endl;
                //std::cout << " ! " << std::flush;
                saddleInt=0;
          }
          maxContrast=float(maxInt-saddleInt);
          maxIntensity=float(maxInt);
 
          std::set<Site,ltstr_p3d<Site> > listePoints=blob->GetListePoints();
          typename std::set<Site,ltstr_p3d<Site> >::iterator
              itPoints=listePoints.begin();
          Site point;
 
           //std::cout << "DEBUG : going through points" << std::endl;
           for (; itPoints!=listePoints.end(); ++itPoints) {
                point=(*itPoints);
                intensity=(*_texdata).intensity(point);
                nbPoint++;
                area=area+1.0;
                meanContrast+=float(intensity - saddleInt);
                meanIntensity+=float(intensity);
           }
           meanContrast/=float(nbPoint);
           meanIntensity/=float(nbPoint);
 
           //std::cout << " + " << std::flush;
           blob->measurements.maxIntensity=maxIntensity;
           blob->measurements.meanIntensity=meanIntensity;
           blob->measurements.maxContrast=maxContrast;
           blob->measurements.meanContrast=meanContrast;
           blob->measurements.area=area;
 
 
          if ( GetOriginalTexture() != NULL ) {
              std::set<Site,ltstr_p3d<Site> >  pixels;
              pixels = blob->GetListePoints();
 
              typename std::set<Site, ltstr_p3d<Site> >::iterator itPix;
              float tvmax = -100.0;
              for ( itPix = pixels.begin() ; itPix != pixels.end() ; itPix++ ) {
                  if ( float ( GetOriginalTexture()->intensity(*itPix) ) > tvmax )
                      tvmax = float ( GetOriginalTexture()->intensity(*itPix) );
              }
              blob->measurements.t = tvmax;
          }
          else  {
              std::cout << "BLOBMEASUREMENT = 0.0" << std::endl;
              blob->measurements.t = 0.0;
          }
 
 
          //std::cout << " - " << std::flush;
          if ( fileStat != 0 )
              fprintf ( fileStat, "%i %.4f %.4f %.4f %.4f %.4f\n", blob->Label(),
                      maxIntensity, meanIntensity, maxContrast, meanContrast, area);
          //std::cout << " | " << std::flush;
      }
      if ( fileStat != 0 )
          fclose(fileStat);
 }
 
 //----------------------------------------------------------------------------------------
 
  // Image specialisation
 
 template <>
 std::vector<float>
 GreyLevelBlobTools<Point3d>::Boundingbox()
 {
 
     std::set<Point3d,ltstr_p3d<Point3d> >  listeP=_blob->GetListePoints();
     std::set<Point3d,ltstr_p3d<Point3d> >::iterator pointIt=listeP.begin();
 
     float x1,x2,y1,y2,z1,z2;
 
     x1=10000.0; y1=10000.0; z1=10000.0;
     x2=-10000.0; y2=-10000.0; z2=-10000.0;
     std::vector<float> vectF(6);
 
     for (; pointIt != listeP.end(); ++pointIt)
     {
          float x=(float) (*pointIt)[0];
          float y=(float) (*pointIt)[1];
          float z=(float) (*pointIt)[2];
          if (x<x1) x1=x;
          if (x>x2) x2=x;
          if (y<y1) y1=y;
          if (y>y2) y2=y;
          if (z<z1) z1=z;
          if (z>z2) z2=z;
     }
 
     vectF[0]=x1; vectF[1]=y1; vectF[2]=z1;
     vectF[3]=x2; vectF[4]=y2; vectF[5]=z2;
 
     return vectF;
 }
 
 template <>
     std::vector<float>
     GreyLevelBlobTools<Point3d >::Barycenter()
     {
       Point3d max= _blob->GetMaximum()->_node;
       std::vector<float> triplet;
       triplet.push_back(max[0]);
       triplet.push_back(max[1]);
       triplet.push_back(max[2]);
       return triplet;
 
     }
 
  // Surface specialisation
 
 template <>
 std::vector<float>
 GreyLevelBlobTools<std::pair<Point3df,uint> >::Boundingbox()
 {
     std::set<std::pair<Point3df,uint> ,ltstr_p3d<std::pair<Point3df,uint> > >  listeP=_blob->GetListePoints();
     std::set<std::pair<Point3df,uint> ,ltstr_p3d<std::pair<Point3df,uint> > >::iterator pointIt=listeP.begin();
 
     float x1,x2,y1,y2,z1,z2;
 
     x1=10000.0; y1=10000.0; z1=10000.0;
     x2=-10000.0; y2=-10000.0; z2=-10000.0;
 
     std::vector<float> vectF(6);
 
     for (; pointIt != listeP.end(); ++pointIt)
     {
          float x=(float) (*pointIt).first[0];
          float y=(float) (*pointIt).first[1];
          float z=(float) (*pointIt).first[2];
          if (x<x1) x1=x;
          if (x>x2) x2=x;
          if (y<y1) y1=y;
          if (y>y2) y2=y;
          if (z<z1) z1=z;
          if (z>z2) z2=z;
     }
 
     vectF[0]=x1; vectF[1]=y1; vectF[2]=z1;
     vectF[3]=x2; vectF[4]=y2; vectF[5]=z2;
 
     return vectF;
 }
 
 
 
template <>
    std::vector<float>
    GreyLevelBlobTools<std::pair<Point3df,uint> >::Barycenter()
{
  Point3df max= _blob->GetMaximum()->_node.first;
  std::vector<float> triplet;
  triplet.push_back(max[0]);
  triplet.push_back(max[1]);
  triplet.push_back(max[2]);
  return triplet;
 
}
}
#endif