convol.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_CONVOL_H
#define AIMS_SIGNALFILTER_CONVOL_H
 
#include <aims/config/aimsalgo_config.h>
#include <cartodata/volume/volume.h>
#include <cartobase/type/limits.h>
#include <float.h>
 

template <class T>
class AimsConvolution
{
public:
 
  AimsConvolution(bool safe=true) : _safe(safe) {
    set_safe_status(safe);
    }
  virtual ~AimsConvolution() { }
 
  carto::VolumeRef< T > doit( carto::rc_ptr<carto::Volume< T > >&,
                              carto::rc_ptr<carto::Volume< T > >& );
 
  void set_safe_status(bool safe=true) {
    if (safe)
      _doit_voxel_method = &AimsConvolution<T>::doit_voxel_safe;
    else
      _doit_voxel_method = &AimsConvolution<T>::doit_voxel_unsafe;
    _safe = safe;
  }
 
protected:
  T doit_voxel( carto::rc_ptr<carto::Volume< T>  >& img,
                carto::rc_ptr<carto::Volume< T > >& kernel,
                int x, int y, int z, int t );
  virtual T doit_voxel_safe( carto::rc_ptr<carto::Volume< T > >& img,
                             carto::rc_ptr<carto::Volume< T > >& kernel,
                             int x, int y, int z, int t);
  virtual T doit_voxel_unsafe( carto::rc_ptr<carto::Volume< T > >& img,
                               carto::rc_ptr<carto::Volume< T > >& kernel,
                               int x, int y, int z, int t);
  T min() const;
  T max() const;
 
  T (AimsConvolution::*_doit_voxel_method) (
    carto::rc_ptr<carto::Volume< T > >& img,
    carto::rc_ptr<carto::Volume< T > >& kernel, int x, int y, int z, int t);
 
  bool _safe;
};

template <class T>
class AimsMaskedConvolution : public AimsConvolution< T >
{
public:
 
  AimsMaskedConvolution( const carto::rc_ptr<carto::Volume< short > >&mask,
                         short sources=0, bool safe=true)
    : AimsConvolution< T >(safe),
      _sources(sources),
      _mask(mask)
  { }
  virtual ~AimsMaskedConvolution() { }
  virtual T doit_voxel_safe( carto::rc_ptr<carto::Volume< T > >& img,
                             carto::rc_ptr<carto::Volume< T > >& kernel,
                             int x, int y, int z, int t);
  virtual T doit_voxel_unsafe( carto::rc_ptr<carto::Volume< T > >& img,
                               carto::rc_ptr<carto::Volume< T > >& kernel,
                               int x, int y, int z, int t);
 
protected:
  short                 _sources;
  const carto::VolumeRef<short>  _mask;
};
 
 
template< class T > inline
carto::VolumeRef< T >
AimsConvolution< T >::doit( carto::rc_ptr<carto::Volume< T > >& img,
                            carto::rc_ptr<carto::Volume< T > >& kernel )
{
  ASSERT( kernel->getSizeT() == 1 );
  if (_safe == false)
    ASSERT( img->getBorders()[0] >= 1 );
  else
    if( img->getBorders()[0] >= 1 )
      set_safe_status(false);
 
  int x, y, z, t;
 
  int dx = img->getSizeX();
  int dy = img->getSizeY();
  int dz = img->getSizeZ();
  int dt = img->getSizeT();
 
  carto::VolumeRef< T > res( dx, dy, dz, dt );
  res.setVoxelSize( img->getVoxelSize() );
 
  for ( t=0; t<dt; t++ )
    for ( z=0; z<dz; z++ )
      for ( y=0; y<dy; y++ )
        for ( x=0; x<dx; x++ )
        {
          T val = doit_voxel(img, kernel, x, y, z, t);
          if ( val > max() )
            val = max();
          else if ( val < min() )
            val = min();
          res( x, y, z, t ) = val;
        }
 
  return res;
}
 
template< class T > inline
T AimsConvolution< T >::doit_voxel( carto::rc_ptr<carto::Volume< T > >& img,
                                    carto::rc_ptr<carto::Volume< T > >& kernel,
                                    int x, int y, int z, int t)
{
  return (this->*_doit_voxel_method)(img, kernel, x, y, z, t);
};
 
template< class T > inline
T AimsConvolution< T >::doit_voxel_safe(
  carto::rc_ptr<carto::Volume< T > >& img,
  carto::rc_ptr<carto::Volume< T > >& kernel, int x, int y, int z, int t)
{
  int i, j, k, idx, idy, idz;
 
  int mx = kernel->getSizeX();
  int my = kernel->getSizeY();
  int mz = kernel->getSizeZ();
 
  int mx2 = mx / 2;
  int my2 = my / 2;
  int mz2 = mz / 2;
 
  int dx = img->getSizeX();
  int dy = img->getSizeY();
  int dz = img->getSizeZ();
 
  T val = (T)0;
 
  for ( k=0; k<mz; k++ )
    for ( j=0; j<my; j++ )
      for ( i=0; i<mx; i++ )
      {
        idx = x - i + mx2;
        idy = y - j + my2;
        idz = z - k + mz2;
 
        if ( idx >= 0 && idy >= 0 && idz >= 0 &&
             idx < dx && idy < dy && idz < dz )
          val += kernel->at( i, j, k ) * img->at( idx, idy, idz, t );
      }
 
  return val;
}
 
template< class T > inline
T AimsConvolution< T >::doit_voxel_unsafe(
  carto::rc_ptr<carto::Volume< T > >& img,
  carto::rc_ptr<carto::Volume< T > >& kernel, int x, int y, int z, int t)
{
  int i, j, k, idx, idy, idz;
 
  int mx = kernel->getSizeX();
  int my = kernel->getSizeY();
  int mz = kernel->getSizeZ();
 
  int mx2 = mx / 2;
  int my2 = my / 2;
  int mz2 = mz / 2;
 
  T val = (T)0;
 
  for ( k=0; k<mz; k++ )
  for ( j=0; j<my; j++ )
  for ( i=0; i<mx; i++ )
  {
    idx = x - i + mx2;
    idy = y - j + my2;
    idz = z - k + mz2;
 
    val += kernel->at( i, j, k ) * img->at( idx, idy, idz, t );
  }
 
  return val;
}
 
template< class T > inline
T AimsMaskedConvolution< T >::doit_voxel_safe(
  carto::rc_ptr<carto::Volume< T > >& img,
  carto::rc_ptr<carto::Volume< T > >& kernel, int x, int y, int z, int t)
{
  if (_mask(x, y, z, t) == _sources)
    return img->at(x, y, z, t);
  return AimsConvolution<T>::doit_voxel_safe(img, kernel, x, y, z, t);
}
 
template< class T > inline
T AimsMaskedConvolution< T >::doit_voxel_unsafe(
  carto::rc_ptr<carto::Volume< T > >& img,
  carto::rc_ptr<carto::Volume< T > >& kernel, int x, int y, int z, int t )
{
  if (_mask(x, y, z, t) == _sources)
    return img->at(x, y, z, t);
  return AimsConvolution<T>::doit_voxel_unsafe(img, kernel, x, y, z, t);
}
 
template <>
inline
int8_t AimsConvolution< int8_t >::min() const
{
  return std::numeric_limits<int8_t>::min();
}
 
 
template <>
inline 
int8_t AimsConvolution< int8_t >::max() const
{
  return std::numeric_limits<int8_t>::max();
}
 
 
template <>
inline
uint8_t AimsConvolution< uint8_t >::min() const
{
  return 0;
}
 
 
template <>
inline 
uint8_t AimsConvolution< uint8_t >::max() const
{
  return std::numeric_limits<uint8_t>::max();
}
 
 
template <>
inline
int16_t AimsConvolution< int16_t >::min() const
{
  return std::numeric_limits<int16_t>::min();
}
 
 
template <>
inline
int16_t AimsConvolution< int16_t >::max() const
{
  return std::numeric_limits<int16_t>::max();
}
 
 
template <>
inline
uint16_t AimsConvolution< uint16_t >::min() const
{
  return 0;
}
 
 
template <>
inline
uint16_t AimsConvolution< uint16_t >::max() const
{
  return std::numeric_limits<uint16_t>::max();
}
 
 
template <>
inline
int32_t AimsConvolution< int32_t >::min() const
{
  return std::numeric_limits<int32_t>::min();
}
 
 
template <>
inline
int32_t AimsConvolution< int32_t >::max() const
{
  return std::numeric_limits<int32_t>::max();
}
 
 
template <>
inline
uint32_t AimsConvolution< uint32_t >::min() const
{
  return 0;
}
 
 
template <>
inline
uint32_t AimsConvolution< uint32_t >::max() const
{
  return std::numeric_limits<uint32_t>::max();
}
 
 
template <>
inline
float AimsConvolution< float >::min() const
{
  return (-FLT_MAX);
}
 
 
template <>
inline
float AimsConvolution< float >::max() const
{
  return FLT_MAX;
}
 
 
template <>
inline
double AimsConvolution< double >::min() const
{
  return DBL_MIN;
}
 
 
template <>
inline
double AimsConvolution< double >::max() const
{
  return DBL_MAX;
}
 
 
// ### remove after everything has been moved to intN_t/uintN_t
#include <limits.h>
#if !defined(__sun__) || !defined(_CHAR_IS_SIGNED)
template<>
inline
char AimsConvolution< char >::min() const
{
  return std::numeric_limits<char>::min();
}
 
 
template <>
inline
char AimsConvolution< char >::max() const
{
  return std::numeric_limits<char>::max();
}
#endif
 
 
template <>
inline
long AimsConvolution< long >::min() const
{
  return LONG_MIN;
}
 
 
template <>
inline
long AimsConvolution< long >::max() const
{
  return LONG_MAX;
}
 
 
template <>
inline
ulong AimsConvolution< ulong >::min() const
{
  return 0;
}
 
 
template <>
inline
ulong AimsConvolution< ulong >::max() const
{
  return ULONG_MAX;
}
 
 
#endif