lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx

   1 // -------------------------------------------------------------------------
   2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
   3 // -------------------------------------------------------------------------
   4
   5 #ifndef __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
   6 #define __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
   7
   8 #include <cstdarg>
   9
  10 // -------------------------------------------------------------------------
  11 template< class _TScalar, unsigned int _VDimension >
  12 void cpExtensions::Algorithms::
  13 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  14 Clear( )
  15 {
  16   this->m_NumberOfSamples = 0;
  17   this->m_Mean.Fill( TScalar( 0 ) );
  18   this->m_Covariance.Fill( TScalar( 0 ) );
  19   this->m_InverseUpdated = false;
  20   this->m_InverseUnbiasedCovariance.Fill( TScalar( 0 ) );
  21 }
  22
  23 // -------------------------------------------------------------------------
  24 template< class _TScalar, unsigned int _VDimension >
  25 template< class _TOtherScalar >
  26 void cpExtensions::Algorithms::
  27 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  28 AddSample( const _TOtherScalar& x1, ... )
  29 {
  30   TVector s;
  31   std::va_list lst;
  32   va_start( lst, x1 );
  33   s[ 0 ] = TScalar( x1 );
  34   for( unsigned int d = 1; d < _VDimension; ++d )
  35     s[ d ] = TScalar( va_arg( lst, double ) );
  36   va_end( lst );
  37   this->_AddSample( s );
  38 }
  39
  40 // -------------------------------------------------------------------------
  41 template< class _TScalar, unsigned int _VDimension >
  42 template< class _TOtherScalar >
  43 void cpExtensions::Algorithms::
  44 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  45 AddSample( const _TOtherScalar* array )
  46 {
  47   TVector s;
  48   for( unsigned d = 0; d < _VDimension; ++d )
  49     s[ d ] = TScalar( array[ d ] );
  50   this->_AddSample( s );
  51 }
  52
  53 // -------------------------------------------------------------------------
  54 template< class _TScalar, unsigned int _VDimension >
  55 template< class _TOtherScalar >
  56 void cpExtensions::Algorithms::
  57 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  58 AddSample( const vnl_vector< _TOtherScalar >& v )
  59 {
  60   unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
  61   TVector s;
  62   for( unsigned d = 0; d < len; ++d )
  63     s[ d ] = TScalar( v[ d ] );
  64   this->_AddSample( s );
  65 }
  66
  67 // -------------------------------------------------------------------------
  68 template< class _TScalar, unsigned int _VDimension >
  69 template< class _TOtherScalar >
  70 void cpExtensions::Algorithms::
  71 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  72 AddSample( const itk::CovariantVector< _TOtherScalar, _VDimension >& v )
  73 {
  74   TVector s;
  75   for( unsigned d = 0; d < _VDimension; ++d )
  76     s[ d ] = TScalar( v[ d ] );
  77   this->_AddSample( s );
  78 }
  79
  80 // -------------------------------------------------------------------------
  81 template< class _TScalar, unsigned int _VDimension >
  82 template< class _TOtherScalar >
  83 void cpExtensions::Algorithms::
  84 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  85 AddSample( const itk::Point< _TOtherScalar, _VDimension >& p )
  86 {
  87   TVector s;
  88   for( unsigned d = 0; d < _VDimension; ++d )
  89     s[ d ] = TScalar( p[ d ] );
  90   this->_AddSample( s );
  91 }
  92
  93 // -------------------------------------------------------------------------
  94 template< class _TScalar, unsigned int _VDimension >
  95 template< class _TOtherScalar >
  96 void cpExtensions::Algorithms::
  97 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
  98 AddSample( const itk::Vector< _TOtherScalar, _VDimension >& v )
  99 {
 100   TVector s;
 101   for( unsigned d = 0; d < _VDimension; ++d )
 102     s[ d ] = TScalar( v[ d ] );
 103   this->_AddSample( s );
 104 }
 105
 106 // -------------------------------------------------------------------------
 107 template< class _TScalar, unsigned int _VDimension >
 108 template< class _TOtherScalar >
 109 _TScalar cpExtensions::Algorithms::
 110 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 111 SquaredMahalanobis( const _TOtherScalar& x1, ... ) const
 112 {
 113   TVector s;
 114   std::va_list lst;
 115   va_start( lst, x1 );
 116   s[ 0 ] = TScalar( x1 );
 117   for( unsigned int d = 1; d < _VDimension; ++d )
 118     s[ d ] = TScalar( va_arg( lst, double ) );
 119   va_end( lst );
 120   return( this->_SquaredMahalanobis( s ) );
 121 }
 122
 123 // -------------------------------------------------------------------------
 124 template< class _TScalar, unsigned int _VDimension >
 125 template< class _TOtherScalar >
 126 _TScalar cpExtensions::Algorithms::
 127 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 128 SquaredMahalanobis( const _TOtherScalar* array ) const
 129 {
 130   TVector s;
 131   for( unsigned d = 0; d < _VDimension; ++d )
 132     s[ d ] = TScalar( array[ d ] );
 133   return( this->_SquaredMahalanobis( s ) );
 134 }
 135
 136 // -------------------------------------------------------------------------
 137 template< class _TScalar, unsigned int _VDimension >
 138 template< class _TOtherScalar >
 139 _TScalar cpExtensions::Algorithms::
 140 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 141 SquaredMahalanobis( const vnl_vector< _TOtherScalar >& v ) const
 142 {
 143   unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
 144   TVector s;
 145   for( unsigned d = 0; d < len; ++d )
 146     s[ d ] = TScalar( v[ d ] );
 147   return( this->_SquaredMahalanobis( s ) );
 148 }
 149
 150 // -------------------------------------------------------------------------
 151 template< class _TScalar, unsigned int _VDimension >
 152 template< class _TOtherScalar >
 153 _TScalar cpExtensions::Algorithms::
 154 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 155 SquaredMahalanobis(
 156   const itk::CovariantVector< _TOtherScalar, _VDimension >& v
 157   ) const
 158 {
 159   TVector s;
 160   for( unsigned d = 0; d < _VDimension; ++d )
 161     s[ d ] = TScalar( v[ d ] );
 162   return( this->_SquaredMahalanobis( s ) );
 163 }
 164
 165 // -------------------------------------------------------------------------
 166 template< class _TScalar, unsigned int _VDimension >
 167 template< class _TOtherScalar >
 168 _TScalar cpExtensions::Algorithms::
 169 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 170 SquaredMahalanobis( const itk::Point< _TOtherScalar, _VDimension >& p ) const
 171 {
 172   TVector s;
 173   for( unsigned d = 0; d < _VDimension; ++d )
 174     s[ d ] = TScalar( p[ d ] );
 175   return( this->_SquaredMahalanobis( s ) );
 176 }
 177
 178 // -------------------------------------------------------------------------
 179 template< class _TScalar, unsigned int _VDimension >
 180 template< class _TOtherScalar >
 181 _TScalar cpExtensions::Algorithms::
 182 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 183 SquaredMahalanobis( const itk::Vector< _TOtherScalar, _VDimension >& v ) const
 184 {
 185   TVector s;
 186   for( unsigned d = 0; d < _VDimension; ++d )
 187     s[ d ] = TScalar( v[ d ] );
 188   return( this->_SquaredMahalanobis( s ) );
 189 }
 190
 191 // -------------------------------------------------------------------------
 192 template< class _TScalar, unsigned int _VDimension >
 193 template< class _TOtherScalar >
 194 _TScalar cpExtensions::Algorithms::
 195 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 196 Density( const _TOtherScalar& x1, ... ) const
 197 {
 198   TVector s;
 199   std::va_list lst;
 200   va_start( lst, x1 );
 201   s[ 0 ] = TScalar( x1 );
 202   for( unsigned int d = 1; d < _VDimension; ++d )
 203     s[ d ] = TScalar( va_arg( lst, double ) );
 204   va_end( lst );
 205   return( this->_Density( s ) );
 206 }
 207
 208 // -------------------------------------------------------------------------
 209 template< class _TScalar, unsigned int _VDimension >
 210 template< class _TOtherScalar >
 211 _TScalar cpExtensions::Algorithms::
 212 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 213 Density( const _TOtherScalar* array ) const
 214 {
 215   TVector s;
 216   for( unsigned d = 0; d < _VDimension; ++d )
 217     s[ d ] = TScalar( array[ d ] );
 218   return( this->_Density( s ) );
 219 }
 220
 221 // -------------------------------------------------------------------------
 222 template< class _TScalar, unsigned int _VDimension >
 223 template< class _TOtherScalar >
 224 _TScalar cpExtensions::Algorithms::
 225 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 226 Density( const vnl_vector< _TOtherScalar >& v ) const
 227 {
 228   unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
 229   TVector s;
 230   for( unsigned d = 0; d < len; ++d )
 231     s[ d ] = TScalar( v[ d ] );
 232   return( this->_Density( s ) );
 233 }
 234
 235 // -------------------------------------------------------------------------
 236 template< class _TScalar, unsigned int _VDimension >
 237 template< class _TOtherScalar >
 238 _TScalar cpExtensions::Algorithms::
 239 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 240 Density(
 241   const itk::CovariantVector< _TOtherScalar, _VDimension >& v
 242   ) const
 243 {
 244   TVector s;
 245   for( unsigned d = 0; d < _VDimension; ++d )
 246     s[ d ] = TScalar( v[ d ] );
 247   return( this->_Density( s ) );
 248 }
 249
 250 // -------------------------------------------------------------------------
 251 template< class _TScalar, unsigned int _VDimension >
 252 template< class _TOtherScalar >
 253 _TScalar cpExtensions::Algorithms::
 254 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 255 Density( const itk::Point< _TOtherScalar, _VDimension >& p ) const
 256 {
 257   TVector s;
 258   for( unsigned d = 0; d < _VDimension; ++d )
 259     s[ d ] = TScalar( p[ d ] );
 260   return( this->_Density( s ) );
 261 }
 262
 263 // -------------------------------------------------------------------------
 264 template< class _TScalar, unsigned int _VDimension >
 265 template< class _TOtherScalar >
 266 _TScalar cpExtensions::Algorithms::
 267 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 268 Density( const itk::Vector< _TOtherScalar, _VDimension >& v ) const
 269 {
 270   TVector s;
 271   for( unsigned d = 0; d < _VDimension; ++d )
 272     s[ d ] = TScalar( v[ d ] );
 273   return( this->_Density( s ) );
 274 }
 275
 276 // -------------------------------------------------------------------------
 277 template< class _TScalar, unsigned int _VDimension >
 278 cpExtensions::Algorithms::
 279 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 280 IterativeGaussianModelEstimator( )
 281   : Superclass( )
 282 {
 283   this->Clear( );
 284 }
 285
 286 // -------------------------------------------------------------------------
 287 template< class _TScalar, unsigned int _VDimension >
 288 cpExtensions::Algorithms::
 289 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 290 ~IterativeGaussianModelEstimator( )
 291 {
 292 }
 293
 294 // -------------------------------------------------------------------------
 295 template< class _TScalar, unsigned int _VDimension >
 296 void cpExtensions::Algorithms::
 297 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 298 _AddSample( const TVector& v ) const
 299 {
 300   this->m_NumberOfSamples += 1;
 301   TScalar a = TScalar( 1 ) / TScalar( this->m_NumberOfSamples );
 302   TScalar b = TScalar( 1 ) - a;
 303   TVector c = v - this->m_Mean;
 304   TMatrix m;
 305   for( unsigned int i = 0; i < _VDimension; ++i )
 306   {
 307     m[ i ][ i ] = c[ i ] * c[ i ];
 308     for( unsigned int j = i + 1; j < _VDimension; ++j )
 309       m[ i ][ j ] = m[ j ][ i ] = c[ i ] * c[ j ];
 310
 311   } // rof
 312   this->m_Covariance = ( this->m_Covariance + ( m * a ) ) * b;
 313   this->m_Mean = ( this->m_Mean * b ) + ( v * a );
 314
 315   TScalar u =
 316     TScalar( this->m_NumberOfSamples ) /
 317     TScalar( this->m_NumberOfSamples - 1 );
 318   this->m_UnbiasedCovariance = this->m_Covariance * u;
 319   this->m_InverseUpdated = false;
 320 }
 321
 322 // -------------------------------------------------------------------------
 323 template< class _TScalar, unsigned int _VDimension >
 324 _TScalar cpExtensions::Algorithms::
 325 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 326 _SquaredMahalanobis( const TVector& v ) const
 327 {
 328   if( !this->m_InverseUpdated )
 329   {
 330     double d =
 331       double(
 332         vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) )
 333         );
 334     if( std::fabs( d ) > double( 0 ) )
 335       this->m_InverseUnbiasedCovariance =
 336         this->m_UnbiasedCovariance.GetInverse( );
 337     else
 338       this->m_InverseUnbiasedCovariance.SetIdentity( );
 339     this->m_InverseUpdated = true;
 340
 341   } // fi
 342   TVector x = v - this->m_Mean;
 343   return( x * ( this->m_InverseUnbiasedCovariance * x ) );
 344 }
 345
 346 // -------------------------------------------------------------------------
 347 template< class _TScalar, unsigned int _VDimension >
 348 _TScalar cpExtensions::Algorithms::
 349 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 350 _Density( const TVector& v ) const
 351 {
 352   static const double _2pik =
 353     std::pow( double( 2 ) * double( vnl_math::pi ), _VDimension );
 354
 355   double s = double( this->_SquaredMahalanobis( v ) ) / double( 2 );
 356   double d =
 357     double(
 358       vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) )
 359       );
 360   return( _TScalar( std::exp( -s ) / std::sqrt( _2pik * d ) ) );
 361 }
 362
 363 #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
 364
 365 // eof - $RCSfile$