iirfilter.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_SIGNALFILTER_IIRFILTER_H
#define AIMS_SIGNALFILTER_IIRFILTER_H

#include <aims/signalfilter/types.h>
#include <cartodata/volume/volume.h>
#include <vector>

namespace aims {

  //==========================================================================
  //
  //   BASE CLASS FOR INFINITE IMPULSE RESPONSE FILTER
  //
  //==========================================================================
  class IIRFilterBase
  {
  public:

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    IIRFilterBase();
    IIRFilterBase( const IIRFilterBase & other );
    virtual ~IIRFilterBase();
    IIRFilterBase & operator=( const IIRFilterBase & other );

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    template <typename T>
    carto::VolumeRef<double> execute( const carto::VolumeRef<T> & in ) const;
    carto::VolumeRef<double> execute( carto::VolumeRef<double> & in, bool inplace = false ) const;
    virtual void filter1d( const carto::VolumeRef<double> in,
                           carto::VolumeRef<double> out, int dir = -1 ) const = 0;

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual const std::vector<bool> & directions() const;
    virtual void setDirections( const std::vector<bool> & dir );
    virtual void setDirections( bool dirx, bool diry, bool dirz, bool dirt );
    virtual void setDirectionX( bool dirx ); 
    virtual void setDirectionY( bool diry ); 
    virtual void setDirectionZ( bool dirz ); 
    virtual void setDirectionT( bool dirt ); 
    virtual const FilterType::Boundary & boundary() const;
    virtual void setBoundary( const FilterType::Boundary & boundary );
    virtual void setVerbose( int verbose = 1 );
    virtual void setQuiet();


  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual void filter( carto::VolumeRef<double> vol ) const;
    virtual void filter1d( carto::VolumeRef<double> vol, int dir = -1,
                           carto::VolumeRef<double> tmp
                              = carto::VolumeRef<double>( (carto::Volume<double>*)0 ) ) const;
    virtual double & at( carto::VolumeRef<double> & vector, long k, int dir ) const;
    virtual const double & at( const carto::VolumeRef<double> & vector, long k, int dir ) const;
    virtual double at( const carto::VolumeRef<double> & vector, long k,
                       int dir, long fullsize, int binf ) const;
    virtual long mirrorCoeff( long i, long size ) const;
    virtual long selectCoeff( long i, long fullsize, int binf ) const;
    virtual long computeK0( double pole, double tolerance = -1. ) const;

    int                   _verbose;
    std::vector<bool>     _dir;
    FilterType::Boundary  _boundary;
    FilterType::CopyType  _copy;
  };

// Yael - April 15, 2016
// I'm commenting out this class for now because the initialization
// strategy is currently false. I guess there must be a right way to do this
// but I don't have the time and do not need this class for now.
// I leave it in case it is of interest someday.
#if 0
  //==========================================================================
  //
  //   CAUSAL INFINITE IMPULSE RESPONSE FILTER
  //
  //==========================================================================
  class CausalIIRFilter: public IIRFilterBase
  {
  public:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    CausalIIRFilter( const std::vector<double> & num_coeffs = std::vector<double>(),
                            const std::vector<double> & den_coeffs = std::vector<double>(),
                            double b0 = 1. );
    CausalIIRFilter( const CausalIIRFilter & other );
    virtual ~CausalIIRFilter();
    CausalIIRFilter & operator=( const CausalIIRFilter & other );

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual void filter1d( const carto::VolumeRef<double> in,
                           carto::VolumeRef<double> out, int dir = -1 ) const;

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    const double & b0() const;
    const std::vector<double> & numCoeffs() const;
    const std::vector<double> & denCoeffs() const;
    void setB0( const double & gain );
    void setNumCoeff( const double & coeff );
    void setNumCoeffs( const std::vector<double> & coeff );
    void setDenCoeff( const double & coeff );
    void setDenCoeffs( const std::vector<double> & coeff );

  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    double recursiveHelper( const carto::VolumeRef<double> in,
                            long k, int dir, long fullsize, int binf ) const;

    double _b0;
    std::vector<double> _num_coeff;
    std::vector<double> _den_coeff;
  };

  //==========================================================================
  //
  //   ALL POLE CAUSAL INFINITE IMPULSE RESPONSE FILTER
  //
  //==========================================================================
  class CausalAllPoleIIRFilter: public IIRFilterBase
  {
  public:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    CausalAllPoleIIRFilter( const std::vector<double> & coeffs = std::vector<double>(),
                            double gain = 1., bool is_pole = false );
    CausalAllPoleIIRFilter( double coeff, double gain = 1., bool is_pole = false );
    CausalAllPoleIIRFilter( const CausalAllPoleIIRFilter & other );
    virtual ~CausalAllPoleIIRFilter();
    CausalAllPoleIIRFilter & operator=( const CausalAllPoleIIRFilter & other );

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual void filter1d( const carto::VolumeRef<double> in,
                           carto::VolumeRef<double> out, int dir = -1 ) const;

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    const double & gain() const;
    const std::vector<double> & poles() const;      
    const std::vector<double> & coeffs() const;     
    void setGain( const double & gain );
    void setPole( const double & pole, bool is_pole = true );
    void setPoles( const std::vector<double> & pole, bool is_pole = true );
    void setCoeff( const double & coeff, bool is_pole = false );
    void setCoeffs( const std::vector<double> & coeff, bool is_pole = false );

  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    void computeK0();

    double _gain;
    std::vector<double> _coeff;
    std::vector<long>   _k0;
    bool   _is_pole;
  };
#endif

  //==========================================================================
  //
  //   SINGLE POLE CAUSAL INFINITE IMPULSE RESPONSE FILTER
  //
  //==========================================================================
  class CausalSinglePoleIIRFilter: public IIRFilterBase
  {
  public:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    CausalSinglePoleIIRFilter( double pole = 0., double gain = 1. );
    CausalSinglePoleIIRFilter( const CausalSinglePoleIIRFilter & other );
    virtual ~CausalSinglePoleIIRFilter();
    CausalSinglePoleIIRFilter & operator=( const CausalSinglePoleIIRFilter & other );

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual void filter1d( const carto::VolumeRef<double> in,
                           carto::VolumeRef<double> out, int dir = -1 ) const;

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    const double & gain() const;
    const double & pole() const;
    void setPole( const double & pole );
    void setGain( const double & gain );

  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    void computeK0();

    double _gain;
    double _pole;
    long   _k0;
  };


// Yael - April 15, 2016
// Same as the CausalFilter class. I think this one is accurate, but it needs
// CausalFilter to work, so it will wait.
#if 0
  //==========================================================================
  //
  //   WARPER FOR ANTICAUSAL FILTERS
  //
  //==========================================================================
  template <typename CausalFilter>
  class AntiCausalWarperIIRFilter: public CausalFilter
  {
  public:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    explicit AntiCausalWarperIIRFilter( const CausalFilter & other = CausalFilter() );
    AntiCausalWarperIIRFilter( const AntiCausalWarperIIRFilter & other );
    virtual ~AntiCausalWarperIIRFilter();
    AntiCausalWarperIIRFilter & operator=( const AntiCausalWarperIIRFilter & other );

  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual double & at( carto::VolumeRef<double> & vector, long k, int dir ) const;
    virtual const double & at( const carto::VolumeRef<double> & vector, long k, int dir ) const;
    virtual double at( const carto::VolumeRef<double> & vector, long k,
                       int dir, long fullsize, int binf ) const;
  };

  //--------------------------------------------------------------------------
  // Typedefs
  //--------------------------------------------------------------------------
  typedef AntiCausalWarperIIRFilter<CausalIIRFilter>           AntiCausalIIRFilter;
  typedef AntiCausalWarperIIRFilter<CausalAllPoleIIRFilter>    AntiCausalAllPoleIIRFilter;
  typedef AntiCausalWarperIIRFilter<CausalSinglePoleIIRFilter> AntiCausalSinglePoleIIRFilter;
#endif

  //==========================================================================
  //
  //   SYMMETRIC ALL-POLE INFINITE IMPULSE RESPONSE FILTER
  //
  //==========================================================================
  class SymAllPoleIIRFilter: public IIRFilterBase
  {
  public:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    SymAllPoleIIRFilter( const double gain = 0.,
                         const std::vector<double> & poles = std::vector<double>() );
    SymAllPoleIIRFilter( const SymAllPoleIIRFilter & other );
    virtual ~SymAllPoleIIRFilter();
    SymAllPoleIIRFilter & operator=( const SymAllPoleIIRFilter & other );

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    virtual void filter1d( const carto::VolumeRef<double> in,
                           carto::VolumeRef<double>       out,
                           int                            dir = -1 ) const;

    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    const double & gain() const;
    const std::vector<double> & poles() const;
    void setPoles( const std::vector<double> & poles );
    void addPole( const double & pole );
    void setGain( const double & gain );
    void setModePreciseGain( bool precise = true );

  protected:
    //------------------------------------------------------------------------
    //------------------------------------------------------------------------
    void computeK0();
    void computePreciseGain();

    std::vector<double>   _poles;
    std::vector<long>     _k0;
    double                _gain;
    double                _precisegain;
    bool                  _precise;
  };

} // namespace aims

#endif // AIMS_SIGNALFILTER_IIRFILTER_H