A.I.M.S algorithms


filteringfunction_linear_d.h
Go to the documentation of this file.
1 /* This software and supporting documentation are distributed by
2  * Institut Federatif de Recherche 49
3  * CEA/NeuroSpin, Batiment 145,
4  * 91191 Gif-sur-Yvette cedex
5  * France
6  *
7  * This software is governed by the CeCILL-B license under
8  * French law and abiding by the rules of distribution of free software.
9  * You can use, modify and/or redistribute the software under the
10  * terms of the CeCILL-B license as circulated by CEA, CNRS
11  * and INRIA at the following URL "http://www.cecill.info".
12  *
13  * As a counterpart to the access to the source code and rights to copy,
14  * modify and redistribute granted by the license, users are provided only
15  * with a limited warranty and the software's author, the holder of the
16  * economic rights, and the successive licensors have only limited
17  * liability.
18  *
19  * In this respect, the user's attention is drawn to the risks associated
20  * with loading, using, modifying and/or developing or reproducing the
21  * software by the user in light of its specific status of free software,
22  * that may mean that it is complicated to manipulate, and that also
23  * therefore means that it is reserved for developers and experienced
24  * professionals having in-depth computer knowledge. Users are therefore
25  * encouraged to load and test the software's suitability as regards their
26  * requirements in conditions enabling the security of their systems and/or
27  * data to be ensured and, more generally, to use and operate it in the
28  * same conditions as regards security.
29  *
30  * The fact that you are presently reading this means that you have had
31  * knowledge of the CeCILL-B license and that you accept its terms.
32  */
33 
34 #ifndef AIMS_SIGNALFILTER_FILTERINGFUNCTION_LINEAR_D_H
35 #define AIMS_SIGNALFILTER_FILTERINGFUNCTION_LINEAR_D_H
36 
37 //--- aims -------------------------------------------------------------------
39 #include <aims/connectivity/structuring_element.h> // aims::StructuringElement
40 #include <aims/io/writer.h>
41 //--- cartodata --------------------------------------------------------------
42 #include <cartodata/volume/volume.h> // carto::VolumeRef
43 //--- cartobase --------------------------------------------------------------
44 #include <cartobase/config/verbose.h> // carto::verbose
45 #include <cartobase/smart/rcptr.h> // carto::rc_ptr
46 //--- std --------------------------------------------------------------------
47 #include <cmath> // std::fabs
48 #include <algorithm> // std::max
49 #include <map>
50 #include <set>
51 #include <string>
52 //----------------------------------------------------------------------------
53 
54 namespace aims {
55 
56 //============================================================================
57 // LINEAR FILTERING FUNCTIONS: FACTORY
58 //============================================================================
59 
60  template <typename T>
62  {
63  static bool initialized = false;
64  if( !initialized )
65  {
66  initialized = true;
67  registerFunction( "gab", GaborFilterFunc<T>() );
68  registerFunction( "gabor", GaborFilterFunc<T>() );
69  }
70  }
71 
72  template <typename T>
73  std::map<std::string,carto::rc_ptr<LinearFilteringFunction<T> > > &
75  {
76  static std::map<std::string,carto::rc_ptr<LinearFilteringFunction<T> > > m;
77  return m;
78  }
79 
80  template <typename T>
82  const std::string & name,
83  const LinearFilteringFunction<T> & func
84  )
85  {
86  init();
87  _map()[ name ] = carto::rc_ptr<LinearFilteringFunction<T> >( func.clone() );
88  }
89 
90  template <typename T>
92  {
93  init();
94  std::set<std::string> s;
95  typename std::map<std::string,carto::rc_ptr<LinearFilteringFunction<T> > >::const_iterator i, e = _map().end();
96  for( i=_map().begin(); i!=e; ++i )
97  s.insert( i->first );
98  return( s );
99  }
100 
101  template <typename T>
103  const std::string & name,
104  carto::Object options
105  )
106  {
107  init();
108  typename std::map<std::string,carto::rc_ptr<LinearFilteringFunction<T> > >::const_iterator i;
109  i = _map().find( name );
110  if( i == _map().end() )
111  return( 0 );
112  LinearFilteringFunction<T> * new_func = i->second->clone();
113  new_func->setOptions( options );
114  return new_func;
115  }
116 
117 //============================================================================
118 // LINEAR FILTERING FUNCTIONS: DERIVED CLASSES
119 //============================================================================
120  //--------------------------------------------------------------------------
121  // GaborFilterFunc
122  //--------------------------------------------------------------------------
123 
124  template <typename T>
125  GaborFilterFunc<T>::GaborFilterFunc( const GaborFilterFunc<T> & other ):
126  _sigma(other._sigma),
127  _theta(other._theta),
128  _lambda(other._lambda),
129  _psi(other._psi),
130  _gamma(other._gamma),
131  _real(other._real),
132  _nstd(other._nstd),
133  _init(other._init),
134  _amplitude(other._amplitude),
135  _voxelsize(other._voxelsize),
136  _kernel(new carto::Volume<double>( *(other._kernel) ))
137  {}
138 
139  template <typename T>
140  GaborFilterFunc<T>::~GaborFilterFunc()
141  {}
142 
143  template <typename T>
144  GaborFilterFunc<T> & GaborFilterFunc<T>::operator=( const GaborFilterFunc<T> & other )
145  {
146  if(this != &other) {
147  _init = other._init ;
148  _sigma = other._sigma;
149  _theta = other._theta;
150  _lambda = other._lambda;
151  _psi = other._psi;
152  _gamma = other._gamma;
153  _real = other._real;
154  _nstd = other._nstd;
155  _voxelsize = other._voxelsize;
156  _amplitude = other._amplitude;
157  _kernel = carto::VolumeRef<double>( new carto::Volume<double>( *(other._kernel) ) );
158  }
159  return *this;
160  }
161 
162  template <typename T>
163  void GaborFilterFunc<T>::setOptions( const carto::Object & options )
164  {
165  _init = false;
166  _sigma = 1.0;
167  _theta = 0.0;
168  _lambda = 1.0;
169  _psi = 0.0;
170  _gamma = 1.0;
171  _real = true;
172  _nstd = 2;
173  _voxelsize = std::vector<float>( 4, 1. );
174  _amplitude = std::vector<int>( 6, 0 );
175 
176  if( options )
177  {
178  if( options->hasProperty( "sigma" ) )
179  _sigma = (double)options->getProperty( "sigma" )->getScalar();
180  if( options->hasProperty( "theta" ) )
181  _theta = (double)options->getProperty( "theta" )->getScalar();
182  if( options->hasProperty( "lambda" ) )
183  _lambda = (double)options->getProperty( "lambda" )->getScalar();
184  if( options->hasProperty( "psi" ) )
185  _psi = (double)options->getProperty( "psi" )->getScalar();
186  if( options->hasProperty( "gamma" ) )
187  _gamma = (double)options->getProperty( "gamma" )->getScalar();
188  if( options->hasProperty( "real" ) )
189  _real = (bool)options->getProperty( "real" )->getScalar();
190  if( options->hasProperty( "nstd" ) )
191  _nstd = (int)options->getProperty( "nstd" )->getScalar();
192  if( options->hasProperty( "voxel_size" ) )
193  options->getProperty( "voxel_size", _voxelsize );
194  }
195 
196  if( carto::verbose )
197  {
198  std::cout << "Gabor Filter created with parameters:" << std::endl
199  << "- sigma: " << _sigma << " mm" << std::endl
200  << "- theta: " << _theta << "°" << std::endl
201  << "- lambda: " << _lambda << " mm" << std::endl
202  << "- psi: " << _psi << "°" << std::endl
203  << "- gamma: " << _gamma << std::endl
204  << "- return " << ( _real ? "real part" : "imaginary part" ) << std::endl
205  << "- nstd: " << _nstd << std::endl
206  << std::endl;
207  }
208 
209  setKernel();
210  }
211 
212  template <typename T>
213  void GaborFilterFunc<T>::updateOptions( const carto::Object & options )
214  {
215  bool changed = false;
216  double oldd;
217  bool oldb;
218  int oldi;
219  std::vector<float> oldv;
220 
221  if( options )
222  {
223  if( options->hasProperty( "sigma" ) ) {
224  oldd = _sigma;
225  _sigma = (double)options->getProperty( "sigma" )->getScalar();
226  if( oldd != _sigma )
227  changed = true;
228  }
229  if( options->hasProperty( "theta" ) ) {
230  oldd = _theta;
231  _theta = (double)options->getProperty( "theta" )->getScalar();
232  if( oldd != _theta )
233  changed = true;
234  }
235  if( options->hasProperty( "lambda" ) ) {
236  oldd = _lambda;
237  _lambda = (double)options->getProperty( "lambda" )->getScalar();
238  if( oldd != _lambda )
239  changed = true;
240  }
241  if( options->hasProperty( "psi" ) ) {
242  oldd = _psi;
243  _psi = (double)options->getProperty( "psi" )->getScalar();
244  if( oldd != _psi )
245  changed = true;
246  }
247  if( options->hasProperty( "gamma" ) ) {
248  oldd = _gamma;
249  _gamma = (double)options->getProperty( "gamma" )->getScalar();
250  if( oldd != _gamma )
251  changed = true;
252  }
253  if( options->hasProperty( "real" ) ) {
254  oldb = _real;
255  _real = (bool)options->getProperty( "real" )->getScalar();
256  if( oldb != _real )
257  changed = true;
258  }
259  if( options->hasProperty( "nstd" ) ) {
260  oldi = _nstd;
261  _nstd = (int)options->getProperty( "nstd" )->getScalar();
262  if( oldi != _nstd )
263  changed = true;
264  }
265  if( options->hasProperty( "voxel_size" ) ) {
266  oldv = _voxelsize;
267  options->getProperty( "voxel_size", _voxelsize );
268  if( oldv != _voxelsize )
269  changed = true;
270  }
271  }
272 
273  if( changed )
274  {
275  if( carto::verbose )
276  std::cout << "Gabor Filter updated with parameters:" << std::endl
277  << "- sigma: " << _sigma << " mm" << std::endl
278  << "- theta: " << _theta << "°" << std::endl
279  << "- lambda: " << _lambda << " mm" << std::endl
280  << "- psi: " << _psi << "°" << std::endl
281  << "- gamma: " << _gamma << std::endl
282  << "- return " << ( _real ? "real part" : "imaginary part" ) << std::endl
283  << "- nstd: " << _nstd << std::endl
284  << std::endl;
285  setKernel();
286  }
287  }
288 
289  template <typename T>
290  const std::vector<int> & GaborFilterFunc<T>::getAmplitude() const
291  {
292  return _amplitude;
293  }
294 
295  template <typename T> inline
296  void GaborFilterFunc<T>::setKernel()
297  {
298  if( carto::verbose )
299  std::cout << "setKernel" << std::endl
300  << "- voxel size: "
301  << _voxelsize[0] << ", "
302  << _voxelsize[1] << std::endl;
303 
304  double pi = 3.1415926535897;
305  std::vector<float> famplitude(2,0.);
306  famplitude[0] = std::max(
307  std::fabs( _nstd * _sigma * std::cos(_theta*pi/180. ) ),
308  std::fabs( _nstd * _sigma/_gamma * std::sin(_theta*pi/180. ) )
309  );
310  famplitude[1] = std::max(
311  std::fabs( _nstd * _sigma * std::sin(_theta*pi/180. ) ),
312  std::fabs( _nstd * _sigma/_gamma * std::cos(_theta*pi/180. ) )
313  );
314  famplitude[0] /= _voxelsize[0];
315  famplitude[1] /= _voxelsize[1];
316  _amplitude[0] = (int) famplitude[0];
317  _amplitude[1] = (int) famplitude[0];
318  _amplitude[2] = (int) famplitude[1];
319  _amplitude[3] = (int) famplitude[1];
320  _kernel = carto::VolumeRef<double>( _amplitude[0] + _amplitude[1] + 1,
321  _amplitude[2] + _amplitude[3] + 1 );
322  if( carto::verbose )
323  std::cout << "- amplitude: "
324  << _amplitude[0] << ", "
325  << _amplitude[1] << ", "
326  << _amplitude[2] << ", "
327  << _amplitude[3] << ", "
328  << _amplitude[4] << ", "
329  << _amplitude[5] << std::endl;
330 
331 
332  double offsetx = (double)(_kernel->getSizeX())/2;
333  double offsety = (double)(_kernel->getSizeY())/2;
334  double xp, yp;
335  for( int32_t y=0; y<_kernel->getSizeY(); ++y )
336  for( int32_t x=0; x<_kernel->getSizeX(); ++x )
337  {
338  xp = ( (double)x-offsetx ) * _voxelsize[0] * std::cos( _theta*pi/180. ) +
339  ( (double)y-offsety ) * _voxelsize[1] * std::sin( _theta*pi/180. );
340  yp = ( -(double)x+offsetx ) * _voxelsize[0] * std::sin( _theta*pi/180. ) +
341  ( (double)y-offsety ) * _voxelsize[1] * std::cos( _theta*pi/180. );
342  _kernel(x,y) = std::exp( ( std::pow(xp,2) +
343  std::pow(_gamma,2) *
344  std::pow(yp,2) ) /
345  ( -2. * std::pow(_sigma,2) ) );
346  if( _real )
347  _kernel(x,y) *= std::cos( 2. * pi * xp / _lambda + _psi * pi / 180. );
348  else
349  _kernel(x,y) *= std::sin( 2. * pi * xp / _lambda + _psi * pi / 180. );
350  }
351 
352  _init = true;
353  }
354 
355  template <typename T> inline
356  T GaborFilterFunc<T>::execute( const carto::VolumeRef<T> & volume ) const
357  {
358  ASSERT( volume->getSizeZ() == 1 );
359  ASSERT( volume->getSizeT() == 1 );
360  double result = 0;
361 
362  for( int32_t y=0; y<volume->getSizeY(); ++y )
363  for( int32_t x=0; x<volume->getSizeX(); ++x )
364  result += _kernel(x,y) * volume(x,y);
365 
366  return (T) result;
367  }
368 
369 } // namespace aims
370 
371 #endif
int getSizeY() const
virtual bool getProperty(const std::string &key, Object &value) const =0
float max(float x, float y)
Definition: thickness.h:97
static std::map< std::string, carto::rc_ptr< LinearFilteringFunction< T > > > & _map()
virtual void setOptions(const carto::Object &options)
Set the parameters of the filters If a parameter value is not set in the options object, a default value must be assigned.
int getSizeZ() const
virtual bool hasProperty(const std::string &key) const =0
virtual LinearFilteringFunction< T > * clone() const =0
clone
int verbose
int getSizeT() const
static LinearFilteringFunction< T > * create(const std::string &name, carto::Object options=carto::none())
static void registerFunction(const std::string &name, const LinearFilteringFunction< T > &func)
Base class for linear filtering functions.
int getSizeX() const
#define ASSERT(EX)
static std::set< std::string > names()