filteringimagealgorithm.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 AIMSALGO_SIGNALFILTER_FILTERINGIMAGEALGORITHM_H
#define AIMSALGO_SIGNALFILTER_FILTERINGIMAGEALGORITHM_H
 
//--- aims -------------------------------------------------------------------
#include <aims/algorithm/imagealgorithm.h>
#include <aims/utility/progress.h>
#include <aims/connectivity/structuring_element.h>
//--- cartodata --------------------------------------------------------------
#include <cartodata/volume/volume.h>
//--- cartobase --------------------------------------------------------------
#include <cartobase/type/string_conversion.h>
#include <cartobase/type/types.h>
#include <cartobase/type/datatypetraits.h>
#include <cartobase/type/datatypeinfo.h>
#include <cartobase/object/object.h>
//--- std --------------------------------------------------------------------
#include <iostream>
#include <assert.h>
//----------------------------------------------------------------------------
 
namespace aims {
  template <typename T> class FilteringImageAlgorithmInterface;
  template <typename T> class LinearFilteringImageAlgorithm;
  template <typename T> class LinearFilteringFunction;
  template <typename T> class ElementFilteringImageAlgorithm;
  template <typename T> class ElementFilteringFunction;
  namespace singlechannel {
    template <typename T> class LinearFilteringImageAlgorithm;
    template <typename T> class ElementFilteringImageAlgorithm;
  }
 
  //==========================================================================
  //   LINEAR FILTERING IMAGE ALGORITHM: MULTI CHANNELS
  //==========================================================================
  template <typename T>
  class LinearFilteringImageAlgorithm: public ImageAlgorithm<T>
  {
    public:
      //----------------------------------------------------------------------
      //   typedef
      //----------------------------------------------------------------------
      typedef T VoxelType;
      typedef typename carto::DataTypeTraits<T>::ChannelType ChannelType;
      typedef singlechannel::LinearFilteringImageAlgorithm<ChannelType> SingleChannelImageAlgorithmType;
 
      //----------------------------------------------------------------------
      //   image algorithm interface: implementation
      //----------------------------------------------------------------------
      virtual void setOptions( const carto::Object & options )
      {
        _lalgo->setOptions( options );
      }
 
      virtual void updateOptions( const carto::Object & options )
      {
        _lalgo->updateOptions( options );
      }
 
      virtual LinearFilteringImageAlgorithm<T> *clone() const
      {
        return new LinearFilteringImageAlgorithm<T>(*this);
      }
 
      //----------------------------------------------------------------------
      //   constructor (default+copy) & destructor
      //----------------------------------------------------------------------
    protected:
      LinearFilteringImageAlgorithm();
      LinearFilteringImageAlgorithm( const LinearFilteringFunction<ChannelType> & f ):
        ImageAlgorithm<T>(),
        _lalgo( SingleChannelImageAlgorithmType( f ).clone() )
      {
        ImageAlgorithm<T>::_algo.reset( _lalgo );
      }
      LinearFilteringImageAlgorithm( const LinearFilteringImageAlgorithm<T> & other ):
        ImageAlgorithm<T>(),
        _lalgo( other._lalgo->clone() )
      {
        ImageAlgorithm<T>::_algo.reset( _lalgo );
      }
      LinearFilteringImageAlgorithm & operator= ( LinearFilteringImageAlgorithm & other )
      {
        _lalgo = other._lalgo->clone();
        ImageAlgorithm<T>::_algo.reset( _lalgo );
        return *this;
      }
    public:
      virtual ~LinearFilteringImageAlgorithm() {}
 
      //----------------------------------------------------------------------
      //   member
      //----------------------------------------------------------------------
    protected:
      SingleChannelImageAlgorithmType * _lalgo;
 
  };
 
  //==========================================================================
  //   ELEMENT FILTERING IMAGE ALGORITHM: MULTI CHANNELS
  //==========================================================================
  template <typename T>
  class ElementFilteringImageAlgorithm: public ImageAlgorithm<T>
  {
    public:
      //----------------------------------------------------------------------
      //   typedef
      //----------------------------------------------------------------------
      typedef T VoxelType;
      typedef typename carto::DataTypeTraits<T>::ChannelType ChannelType;
      typedef singlechannel::ElementFilteringImageAlgorithm<ChannelType> SingleChannelImageAlgorithmType;
 
      //----------------------------------------------------------------------
      //   image algorithm interface: implementation
      //----------------------------------------------------------------------
      virtual void setOptions( const carto::Object & options )
      {
        _ealgo->setOptions( options );
      }
 
      virtual void updateOptions( const carto::Object & options )
      {
        _ealgo->updateOptions( options );
      }
 
      virtual ElementFilteringImageAlgorithm<T> *clone() const {
        return new ElementFilteringImageAlgorithm<T>(*this);
      }
 
      //----------------------------------------------------------------------
      //   member method
      //----------------------------------------------------------------------
      virtual void setStructuringElement( const StructuringElement & se )
      {
        _ealgo->setStructuringElement( se );
      }
 
      //----------------------------------------------------------------------
      //   constructor (default+copy) & destructor
      //----------------------------------------------------------------------
    protected:
      ElementFilteringImageAlgorithm(
        const ElementFilteringFunction<ChannelType> & f,
        const StructuringElement & se = strel::Cube(1.)
      ):
        ImageAlgorithm<T>(),
        _ealgo( SingleChannelImageAlgorithmType( f, se ).clone() )
      {
        ImageAlgorithm<T>::_algo.reset( _ealgo );
      }
      ElementFilteringImageAlgorithm( const ElementFilteringImageAlgorithm<T> & other ):
        ImageAlgorithm<T>(),
        _ealgo(other._ealgo->clone())
      {
        ImageAlgorithm<T>::_algo.reset( _ealgo );
      }
      ElementFilteringImageAlgorithm & operator= ( ElementFilteringImageAlgorithm & other )
      {
        _ealgo = other._ealgo->clone();
        ImageAlgorithm<T>::_algo.reset( _ealgo );
        return *this;
      }
    public:
      virtual ~ElementFilteringImageAlgorithm() {}
 
      //----------------------------------------------------------------------
      //   member
      //----------------------------------------------------------------------
    protected:
      SingleChannelImageAlgorithmType * _ealgo;
  };
 
  namespace singlechannel {
 
    //========================================================================
    //   LINEAR FILTERING IMAGE ALGORITHM: SINGLE CHANNEL
    //========================================================================
    template <typename T>
    class LinearFilteringImageAlgorithm: public ImageAlgorithmInterface<T>
    {
      public:
        //--------------------------------------------------------------------
        //   image algorithm interface: implementation
        //--------------------------------------------------------------------
        virtual void execute( const carto::VolumeRef<T> & in,
                                    carto::VolumeRef<T> & out ) const;
 
        virtual void setOptions( const carto::Object & options )
        {
          _func->setOptions( options );
        }
 
        virtual void updateOptions( const carto::Object & options )
        {
          _func->updateOptions( options );
        }
 
        virtual ::aims::singlechannel::LinearFilteringImageAlgorithm<T> * clone() const {
          return new ::aims::singlechannel::LinearFilteringImageAlgorithm<T>(*this);
        }
 
        //--------------------------------------------------------------------
        //   constructor (default+copy) & destructor
        //--------------------------------------------------------------------
      protected:
        LinearFilteringImageAlgorithm();
 
      public:
        LinearFilteringImageAlgorithm( const LinearFilteringFunction<T> & f ):
          ImageAlgorithmInterface<T>(),
          _func(f.clone())
        {}
 
        LinearFilteringImageAlgorithm( const ::aims::singlechannel::LinearFilteringImageAlgorithm<T> & other ):
          ImageAlgorithmInterface<T>(),
          _func(other._func->clone())
        {}
 
        LinearFilteringImageAlgorithm & operator= (
          const LinearFilteringImageAlgorithm & other )
        {
          assert( typeid(other) == typeid(*this) );
          if( this != &other )
            _func.reset( other._func->clone() );
          return *this;
        }
 
        virtual ~LinearFilteringImageAlgorithm() {}
 
        //--------------------------------------------------------------------
        //   member
        //--------------------------------------------------------------------
      protected:
        carto::rc_ptr<LinearFilteringFunction<T> > _func;
    };
 
    template <typename T>
    inline
    void LinearFilteringImageAlgorithm<T>::execute(
      const carto::VolumeRef<T> & in,
            carto::VolumeRef<T> & out ) const
    {
      int sx, sy, sz, ex, ey, ez, x, y, z, t;
      typename carto::Volume<T>::Position4Di pos(0,0,0,0), size(1,1,1,1);
 
      // set voxel size
      carto::Object options = carto::Object::value( carto::PropertySet() );
      std::vector<float> vs(4, 1.);
      in.header().getProperty( "voxel_size", vs );
      options->setProperty( "voxel_size", vs );
      _func->updateOptions( options );
 
      // Manage borders / structuring element size
      const std::vector<int> & border = in->getBorders();
      const std::vector<int> & amplitude = _func->getAmplitude();
 
      // When volume has borders, it is possible to use it to process filtering
      sz = (amplitude[4] < border[4]) ? 0 : amplitude[4] - border[4];
      sy = (amplitude[2] < border[2]) ? 0 : amplitude[2] - border[2];
      sx = (amplitude[0] < border[0]) ? 0 : amplitude[0] - border[0];
      ez = (amplitude[5] < border[5]) ? in->getSizeZ() : in->getSizeZ() - amplitude[5] + border[5];
      ey = (amplitude[3] < border[3]) ? in->getSizeY() : in->getSizeY() - amplitude[3] + border[3];
      ex = (amplitude[1] < border[1]) ? in->getSizeX() : in->getSizeX() - amplitude[1] + border[1];
 
      if( ImageAlgorithmInterface<T>::_verbose > 1 ) {
        std::cout << "Filter amplitude (voxels): [ ";
        for(int i = 0; i < 6; ++i)
          std::cout << carto::toString(amplitude[i]) << ( i==5 ? " " : ", ");
        std::cout << "]" << std::endl;
        std::cout << "Processing with borders (voxels): [ ";
        for(int i = 0; i < 8; ++i)
          std::cout << carto::toString(border[i]) << ( i==7 ? " " : ", ");
        std::cout << "]" << std::endl;
        std::cout << "Start: [" << carto::toString(sx) << ", "
                                << carto::toString(sy) << ", "
                                << carto::toString(sz) << "]" << std::endl;
        std::cout << "End: [" << carto::toString(ex) << ", "
                              << carto::toString(ey) << ", "
                              << carto::toString(ez) << "]" << std::endl;
      }
 
      size[0] = amplitude[0] + amplitude[1] + 1;
      size[1] = amplitude[2] + amplitude[3] + 1;
      size[2] = amplitude[4] + amplitude[5] + 1;
      carto::VolumeRef<T> win( new carto::Volume<T>( in, pos, size ) );
 
      aims::Progression progress( static_cast<size_t>(out->getSizeT())
                                  * (ez - sz) * (ey - sy) * (ex - sx) );
      if( ImageAlgorithmInterface<T>::_verbose > 0 )
        std::cout << "Filtering progress: ";
 
      for( t = 0; t < in->getSizeT(); ++t )
        for( z = sz; z < ez; ++z )
          for( y = sy; y < ey; ++y )
            for( x = sx; x < ex; ++x )
            {
              if( this->_verbose > 0 )
                (++progress).print();
              // Set view position in the volume
              pos[0] = x - amplitude[0];
              pos[1] = y - amplitude[2];
              pos[2] = z - amplitude[4];
              pos[3] = t;
              win->setPosInRefVolume( pos );
              (*out)( x, y, z, t ) = _func->execute( win );
            }
 
      if( ImageAlgorithmInterface<T>::_verbose > 0 )
        std::cout << std::endl;
    }
 
 
    //========================================================================
    //   ELEMENT FILTERING IMAGE ALGORITHM: SINGLE CHANNEL
    //========================================================================
    template <typename T>
    class ElementFilteringImageAlgorithm: public ImageAlgorithmInterface<T>
    {
      public:
        //--------------------------------------------------------------------
        //   image algorithm interface: implementation
        //--------------------------------------------------------------------
        virtual void execute( const carto::VolumeRef<T> & in,
                                    carto::VolumeRef<T> & out ) const;
 
        virtual void setOptions( const carto::Object & options )
        {
          _func->setOptions( options );
        }
 
        virtual void updateOptions( const carto::Object & options )
        {
          _func->updateOptions( options );
        }
 
        virtual ::aims::singlechannel::ElementFilteringImageAlgorithm<T> * clone() const {
          return new ::aims::singlechannel::ElementFilteringImageAlgorithm<T>(*this);
        }
 
        //--------------------------------------------------------------------
        //   member method
        //--------------------------------------------------------------------
        virtual void setStructuringElement( const StructuringElement & se )
        {
          _strel.reset( se.clone() );
        }
 
        //--------------------------------------------------------------------
        //   constructor (default+copy) & destructor
        //--------------------------------------------------------------------
      protected:
        ElementFilteringImageAlgorithm();
 
      public:
        ElementFilteringImageAlgorithm( const ElementFilteringFunction<T> & f,
                                        const StructuringElement & se = strel::Cube(1.) ):
          ImageAlgorithmInterface<T>(),
          _func(f.clone()),
          _strel(se.clone())
        {}
 
        ElementFilteringImageAlgorithm( const ::aims::singlechannel::ElementFilteringImageAlgorithm<T> & other ):
          ImageAlgorithmInterface<T>(),
          _func(other._func->clone()),
          _strel(other._strel->clone())
        {}
 
        ::aims::singlechannel::ElementFilteringImageAlgorithm<T> & operator=(
          const ::aims::singlechannel::ElementFilteringImageAlgorithm<T> & other )
        {
          assert( typeid(other) == typeid(*this) );
          if( *this != other )
          {
            _func.reset( other._func->clone() );
            _strel.reset( other._strel->clone() );
          }
          return *this;
        }
 
        virtual ~ElementFilteringImageAlgorithm() {}
 
        //--------------------------------------------------------------------
        //   members
        //--------------------------------------------------------------------
      protected:
        carto::rc_ptr<ElementFilteringFunction<T> > _func;
        carto::rc_ptr<StructuringElement>           _strel;
    };
 
    template <typename T>
    inline
    void ElementFilteringImageAlgorithm<T>::execute(
      const carto::VolumeRef<T> & in,
            carto::VolumeRef<T> & out
    ) const
    {
      int sx, sy, sz, ex, ey, ez, x, y, z, t;
      typename carto::Volume<T>::Position4Di pos(0,0,0), size(1,1,1);
 
      // Manage borders / structuring element size
      const std::vector<int> & border = in->getBorders();
      const std::vector<int> amplitude = _strel->getAmplitude();
 
      // When volume has borders, it is possible to use it to process filtering
      sz = (amplitude[4] < border[4]) ? 0 : amplitude[4] - border[4];
      sy = (amplitude[2] < border[2]) ? 0 : amplitude[2] - border[2];
      sx = (amplitude[0] < border[0]) ? 0 : amplitude[0] - border[0];
      ez = (amplitude[5] < border[5]) ? in->getSizeZ() : in->getSizeZ() - amplitude[5] + border[5];
      ey = (amplitude[3] < border[3]) ? in->getSizeY() : in->getSizeY() - amplitude[3] + border[3];
      ex = (amplitude[1] < border[1]) ? in->getSizeX() : in->getSizeX() - amplitude[1] + border[1];
 
      if( ImageAlgorithmInterface<T>::_verbose > 1 ) {
        std::cout << "Filter amplitude (voxels): [ ";
        for(int i = 0; i < 6; ++i)
          std::cout << carto::toString(amplitude[i]) << ( i==8 ? " " : ", ");
        std::cout << "]" << std::endl;
        std::cout << "Processing with borders (voxels): [ ";
        for(int i = 0; i < 8; ++i)
          std::cout << carto::toString(border[i]) << ( i==8 ? " " : ", ");
        std::cout << "]" << std::endl;
        std::cout << "Start: [" << carto::toString(sx) << ", "
                                << carto::toString(sy) << ", "
                                << carto::toString(sz) << "]" << std::endl;
        std::cout << "End: [" << carto::toString(ex) << ", "
                              << carto::toString(ey) << ", "
                              << carto::toString(ez) << "]" << std::endl;
      }
 
      carto::VolumeRef<T> win( new carto::Volume<T>( in, pos, size ) );
 
      aims::Progression progress( static_cast<size_t>(out->getSizeT())
                                  * (ez - sz) * (ey - sy) *(ex - sx) );
      if( ImageAlgorithmInterface<T>::_verbose > 0 )
        std::cout << "Filtering progress: ";
 
      for( t = 0; t < in->getSizeT(); ++t )
        for( z = sz; z < ez; ++z )
          for( y = sy; y < ey; ++y )
            for( x = sx; x < ex; ++x )
            {
              if( this->_verbose > 0 )
                (++progress).print();
              // Set view position in the volume
              pos[0] = x;
              pos[1] = y;
              pos[2] = z;
              pos[3] = t;
              win->setPosInRefVolume( pos );
              (*out)( x, y, z, t ) = _func->execute( win, _strel );
            }
 
      if( ImageAlgorithmInterface<T>::_verbose > 0 )
        std::cout << std::endl;
    }
 
  } // namespace singlechannel
 
} // namespace aims
 
#endif