coustywatershed.h

Go to the documentation of this file.

/* Copyright (C) 2000-2013 CEA
 *
 * This software and supporting documentation were developed by
 *     bioPICSEL
 *     CEA/DSV/I²BM/MIRCen/LMN, Batiment 61,
 *     18, route du Panorama
 *     92265 Fontenay-aux-Roses
 *     France
 */

#ifndef BIOPROCESSING_WATERSHED_COUSTY
#define BIOPROCESSING_WATERSHED_COUSTY

//--- bioprocessing ----------------------------------------------------------
#include <bioprocessing/watershed/coustyflowmap.h>    // bio::CoustyflowMapRef
#include <bioprocessing/watershed/coustystream.h>       // bio::CoustyStremRef
#include <bioprocessing/graph/volumegraph.h>            // bio::VolumeGraphRef
//--- aims -------------------------------------------------------------------
#include <aims/vector/vector.h>                                      // Point*
//--- cartodata --------------------------------------------------------------
#include <cartodata/volume/volume.h>                       // carto::VolumeRef
//--- cartobase --------------------------------------------------------------
#include <cartobase/config/verbose.h>                        // carto::verbose
//--- std --------------------------------------------------------------------
#include <iostream>                                    // std::cout, std::endl
//----------------------------------------------------------------------------

namespace bio {
  //==========================================================================
  //   COUSTY'S WATERSHED (Cousty et al, IEEE PAMI 2007)
  //==========================================================================
  template <typename T, typename L>
  class CoustyWatershed
  {
    //--- typedef ------------------------------------------------------------
  public:
    typedef Point4dl                      Point;        
    typedef VolumeGraphRef<T,Point>       Graph;        
    typedef typename Graph::Vertex        Vertex;       
    typedef typename Graph::Edge          Edge;         
    typedef CoustyStreamRef<Vertex>       Stream;       
    typedef CoustyFlowMapRef<Vertex,L>    FlowMap;      

    //--- constructor --------------------------------------------------------
  public:
    CoustyWatershed(): _g(), _s(), _psi(), _verbose(carto::verbose), _3d(true) {}
    ~CoustyWatershed() {}

    //--- properties ---------------------------------------------------------
  public:
    void setVerbose( int verbose = 1 ) { _verbose = verbose; }
    void setQuiet() { setVerbose(0); }
    void set3D() { _3d = true; }
    void set2D() { _3d = false; }

    //--- implementation -----------------------------------------------------
  public:
    carto::VolumeRef<L> execute( carto::VolumeRef<T> in );
    void execute( carto::VolumeRef<T> in, carto::VolumeRef<L> out );
    carto::VolumeRef<L> getMinima() { return _psi.getMinima(); }

    //--- helpers ------------------------------------------------------------
  protected:
    void computeWatershed();
    L computeStream( Vertex x );
    bool findEdge( Vertex y, Vertex & z );
    void setLabel( L lab );

    //--- members ------------------------------------------------------------
  protected:
    Graph    _g;
    Stream   _s;
    FlowMap  _psi;
    int      _verbose;
    bool     _3d;
  };

  //==========================================================================
  //   DEFINITION
  //==========================================================================
  //
  template <typename T, typename L>
  carto::VolumeRef<L> CoustyWatershed<T,L>::execute( carto::VolumeRef<T> in )
  {
    carto::VolumeRef<L> out( in.getSizeX(), in.getSizeY(),
                             in.getSizeZ(), in.getSizeT() );
    out->copyHeaderFrom( in.header() );
    execute( in, out );
    return out;
  }

  template <typename T, typename L>
  void CoustyWatershed<T,L>::execute( carto::VolumeRef<T> in,
                                      carto::VolumeRef<L> out )
  {
    if( _3d )
      _g = Graph( in, aims::strel::Connectivity6XYZ() );
    else
    _g = Graph( in, aims::strel::Connectivity4XY() );
    _s = Stream();
    _psi.setVolume( out );
    _psi.init( _g.beginVertex(), _g.endVertex() );
    computeWatershed();
  }

  template <typename T, typename L>
  void CoustyWatershed<T,L>::computeWatershed()
  {
    Vertex x;
    L nb_labs = 0;
    L lab;
    typename FlowMap::iterator i;

    int j=0;
    long size = _psi.size();
    if( _verbose )
      std::cout << std::endl;
    for( i = _psi.begin() ; i != _psi.end(); ++i )
    {
      if( _verbose && ( ++j % 10000 == 0 ) )
        std::cout << "\rVertex: " << j << "/" << size << std::flush;
      x = i->first;
      lab = i->second;
      if( lab == _psi.labelNone() )
      {
        lab = computeStream( x );
        if( lab == _psi.labelMinus() )
          setLabel( ++nb_labs );
        else {
          setLabel( lab );
        }
      }
    }
    if( _verbose )
      std::cout << std::endl << "Labels created: " << nb_labs << std::endl;
  }

  template <typename T, typename L>
  void CoustyWatershed<T,L>::setLabel( L lab )
  {
    typename Stream::iterator i;
    for( i = _s.begin(); i != _s.end(); ++i )
    {
      _psi.setLabel( *i, lab );
    }
  }

  template <typename T, typename L>
  L CoustyWatershed<T,L>::computeStream( Vertex x )
  {
    _s.reset( x );
    Stream s2( x );
    bool breadth_first;
    Vertex y, z;
    L lab;

    while( !s2.empty() )
    {
      y = *(s2.begin());
      s2.erase( s2.begin() );
      breadth_first = true;
      while( breadth_first && findEdge( y, z ) )
      {
        lab = _psi.label(z);
        if( lab != _psi.labelNone() )
          return lab;
        else if( _g.Fm( z ) < _g.Fm( y ) )
        {
          _s.insert( z );
          s2.reset( z );
          breadth_first = false;
        }
        else
        {
          _s.insert( z );
          s2.insert( z );
        }
      }
    }
    _psi.setMinimum(y);
    return _psi.labelMinus();
  }

  template <typename T, typename L>
  bool CoustyWatershed<T,L>::findEdge( Vertex y, Vertex & z )
  {
    typename Graph::edge_iterator e;
    for( e = _g.beginEdge(y); e != _g.endEdge(y); ++e )
    {
      z = e->other( y );
      if( !_s.count( z ) )
      {
        if( _g.weight( *e ) == _g.Fm( y ) )
          return true;
      }
    }
    z = Vertex::none();
    return false;
  }

} // namesapce bio

#endif // BIOPROCESSING_WATERSHED_COUSTY