extremafinder.h

Go to the documentation of this file.

#ifndef BIOPROCESSING_UTILITY_EXTREMAFINDER_H
#define BIOPROCESSING_UTILITY_EXTREMAFINDER_H
 
//============================================================================
// EXTREMA FINDER
//============================================================================
 
//--- aims -------------------------------------------------------------------
#include <aims/connectivity/structuring_element.h>       // StructuringElement
#include <aims/data/cartodatavolume.h>                             // AimsData
#include <aims/bucket/bucketMap.h>                                // BucketMap
#include <aims/utility/progress.h>
//--- cartodata --------------------------------------------------------------
#include <cartodata/volume/volume.h>                              // VolumeRef
//--- cartobase --------------------------------------------------------------
#include <cartobase/config/verbose.h>                               // verbose
#include <cartobase/smart/rcptr.h>                                   // rc_ptr
//--- std --------------------------------------------------------------------
#include <iostream>
//----------------------------------------------------------------------------
 
namespace bio {
 
  class ExtremaFinder
  {
  public:
    ExtremaFinder():
      _strel( new aims::strel::Sphere(1.0) ),
      _connex( new aims::strel::Connectivity26XYZ() ),
      _modemin(true),
      _modestrict(false),
      _value(1),
      _modecc(false)
    {}
    virtual ~ExtremaFinder() {}
 
    //--- parameters
    void setStructuringElement( const aims::StructuringElement & se )
    {
      _strel.reset( se.clone() );
    }
    void setModeMin() { _modemin = true; _modestrict = false; }
    void setModeMax() { _modemin = false; _modestrict = false; }
    void setModeStrictMin() { _modemin = true; _modestrict = true; }
    void setModeStrictMax() { _modemin = false; _modestrict = true; }
    void setModeValueConst( double value ) { _value = value; _modecc = false; }
    void setModeValueConnectComp( const aims::strel::Connectivity & c )
    {
      _connex.reset( c.clone() );
      _modecc = true;
    }
 
    //--- execute
    template <typename IN, typename OUT>
    void computeExtrema( carto::VolumeRef<IN> in,
                         carto::VolumeRef<OUT> out );
    template <typename IN, typename OUT, typename M>
    void computeExtrema( carto::VolumeRef<IN> in,
                         carto::VolumeRef<OUT> out,
                         carto::VolumeRef<M> mask );
    template <typename IN, typename OUT>
    void computeExtrema( carto::VolumeRef<IN> in,
                         aims::BucketMap<OUT> & out );
    template <typename IN, typename OUT, typename M>
    void computeExtrema( carto::VolumeRef<IN> in,
                         aims::BucketMap<OUT> & out,
                         carto::VolumeRef<M> mask );
 
    //--- utility methods
    template <typename IN, typename OUT>
    carto::VolumeRef<OUT> computeExtremaVolume( carto::VolumeRef<IN> in )
    {
      carto::VolumeRef<OUT> out( in.getSizeX(),
                                 in.getSizeY(),
                                 in.getSizeZ(),
                                 in.getSizeT() );
      out->copyHeaderFrom( in.header() );
      out.fill(0);
      computeExtrema( in, out );
      return out;
    }
 
    template <typename IN, typename OUT, typename M >
    carto::VolumeRef<OUT> computeExtremaVolume(
      carto::VolumeRef<IN> in,
      carto::VolumeRef<M> mask )
    {
      carto::VolumeRef<OUT> out( in.getSizeX(),
                                 in.getSizeY(),
                                 in.getSizeZ(),
                                 in.getSizeT() );
      out->copyHeaderFrom( in.header() );
      out.fill(0);
      computeExtrema( in, out, mask );
      return out;
    }
 
    template <typename IN, typename OUT>
    aims::BucketMap<OUT> computeExtremaBucket( carto::VolumeRef<IN> in )
    {
      aims::BucketMap<OUT> out;
      std::vector<float> vs(4,1.);
      if( in.header().hasProperty( "voxel_size" ) )
        in.header().getProperty( "voxel_size", vs );
      out.setSizeXYZT( vs[0], vs[1], vs[2], vs[3] );
      computeExtrema( in, out );
      return out;
    }
 
    template <typename IN, typename OUT, typename M>
    aims::BucketMap<OUT> computeExtremaBucket(
      carto::VolumeRef<IN> in,
      carto::VolumeRef<M> mask
      = carto::VolumeRef<M>( (carto::Volume<M>*)0 ) )
    {
      aims::BucketMap<OUT> out;
      std::vector<float> vs(4,1.);
      if( in.header().hasProperty( "voxel_size" ) )
        in.header().getProperty( "voxel_size", vs );
      out.setSizeXYZT( vs[0], vs[1], vs[2], vs[3] );
      computeExtrema( in, out, mask );
      return out;
    }
 
  private:
    template <typename IN>
    bool isExtremum( carto::VolumeRef<IN> in, int x, int y, int z, int t );
    carto::rc_ptr<aims::StructuringElement> _strel;
    carto::rc_ptr<aims::strel::Connectivity> _connex;
    bool   _modemin;
    bool   _modestrict;
    bool   _modecc;
    double _value;
  };
 
  template <typename IN>
  bool ExtremaFinder::isExtremum( carto::VolumeRef<IN> in,
                                  int x, int y, int z, int t )
  {
    aims::StructuringElement::const_iterator p;
 
    for( p = _strel->begin(); p != _strel->end(); ++p )
    {
      if( _modemin && _modestrict )
      {
        if( in( x, y, z, t ) >= in( x + (*p)[0], y + (*p)[1], z + (*p)[2], t ) )
          return false;
      }
      else if( _modemin && !_modestrict )
      {
        if( in( x, y, z, t ) > in( x + (*p)[0], y + (*p)[1], z + (*p)[2], t ) )
          return false;
      }
      else if( !_modemin && _modestrict )
      {
        if( in( x, y, z, t ) <= in( x + (*p)[0], y + (*p)[1], z + (*p)[2], t ) )
          return false;
      }
      else if( !_modemin && !_modestrict )
      {
        if( in( x, y, z, t ) < in( x + (*p)[0], y + (*p)[1], z + (*p)[2], t ) )
          return false;
      }
    }
 
    return true;
  }
 
  template <typename IN, typename OUT>
  void ExtremaFinder::computeExtrema( carto::VolumeRef<IN> in,
                                      carto::VolumeRef<OUT> out )
  {
    int x, y, z, t;
    for( t=0; t<in.getSizeT(); ++t )
      for( z=0; z<in.getSizeZ(); ++z )
        for( y=0; y<in.getSizeY(); ++y )
          for( x=0; x<in.getSizeX(); ++x )
            if( isExtremum( in, x, y, z, t ) )
              out( x, y, z, t ) = _value;
 
    if( _modecc )
      std::cout << "ExtremaFinder:: connected components mode not implemented yet." << std::endl;
  }
 
  template <typename IN, typename OUT, typename M>
  void ExtremaFinder::computeExtrema( carto::VolumeRef<IN> in,
                                      carto::VolumeRef<OUT> out,
                                      carto::VolumeRef<M> mask )
  {
    int x, y, z, t;
    for( t=0; t<in.getSizeT(); ++t )
      for( z=0; z<in.getSizeZ(); ++z )
        for( y=0; y<in.getSizeY(); ++y )
          for( x=0; x<in.getSizeX(); ++x )
            if( mask( x, y, z, t ) )
              if( isExtremum( in, x, y, z, t ) )
                out( x, y, z, t ) = _value;
 
    if( _modecc )
      std::cout << "ExtremaFinder:: connected components mode not implemented yet." << std::endl;
  }
 
  template <typename IN, typename OUT>
  void ExtremaFinder::computeExtrema( carto::VolumeRef<IN> in,
                                      aims::BucketMap<OUT> & out )
  {
    int x, y, z, t;
 
    double lmax = in.getSizeT();
    lmax *= in.getSizeZ();
    lmax *= in.getSizeY();
    double l=0;
 
    for( t=0; t<in.getSizeT(); ++t )
      for( z=0; z<in.getSizeZ(); ++z )
        for( y=0; y<in.getSizeY(); ++y )
          for( x=0; x<in.getSizeX(); ++x )
            if( isExtremum( in, x, y, z, t ) )
              out.insert( Point3d(x,y,z), (OUT)_value );
 
    if( _modecc )
      std::cout << "ExtremaFinder:: connected components mode not implemented yet." << std::endl;
  }
 
  template <typename IN, typename OUT, typename M>
  void ExtremaFinder::computeExtrema( carto::VolumeRef<IN> in,
                                      aims::BucketMap<OUT> & out,
                                      carto::VolumeRef<M> mask )
  {
    int x, y, z, t;
 
    double lmax = in.getSizeT();
    lmax *= in.getSizeZ();
    lmax *= in.getSizeY();
    double l=0;
 
    for( t=0; t<in.getSizeT(); ++t )
      for( z=0; z<in.getSizeZ(); ++z )
        for( y=0; y<in.getSizeY(); ++y )
          for( x=0; x<in.getSizeX(); ++x )
            if( mask( x, y, z, t ) )
              if( isExtremum( in, x, y, z, t ) )
                out.insert( Point3d(x,y,z), (OUT)_value );
 
    if( _modecc )
      std::cout << "ExtremaFinder:: connected components mode not implemented yet." << std::endl;
  }
 
} // namespace bio
 
#endif