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 cpExtensions::Algorithms::
 194 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 195 IterativeGaussianModelEstimator( )
 196   : Superclass( )
 197 {
 198   this->Clear( );
 199 }
 200
 201 // -------------------------------------------------------------------------
 202 template< class _TScalar, unsigned int _VDimension >
 203 cpExtensions::Algorithms::
 204 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 205 ~IterativeGaussianModelEstimator( )
 206 {
 207 }
 208
 209 // -------------------------------------------------------------------------
 210 template< class _TScalar, unsigned int _VDimension >
 211 void cpExtensions::Algorithms::
 212 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 213 _AddSample( const TVector& v ) const
 214 {
 215   this->m_NumberOfSamples += 1;
 216   TScalar a = TScalar( 1 ) / TScalar( this->m_NumberOfSamples );
 217   TScalar b = TScalar( 1 ) - a;
 218   TVector c = v - this->m_Mean;
 219   TMatrix m;
 220   for( unsigned int i = 0; i < _VDimension; ++i )
 221   {
 222     m[ i ][ i ] = c[ i ] * c[ i ];
 223     for( unsigned int j = i + 1; j < _VDimension; ++j )
 224       m[ i ][ j ] = m[ j ][ i ] = c[ i ] * c[ j ];
 225
 226   } // rof
 227   this->m_Covariance = ( this->m_Covariance + ( m * a ) ) * b;
 228   this->m_Mean = ( this->m_Mean * b ) + ( v * a );
 229
 230   TScalar u =
 231     TScalar( this->m_NumberOfSamples ) /
 232     TScalar( this->m_NumberOfSamples - 1 );
 233   this->m_UnbiasedCovariance = this->m_Covariance * u;
 234   this->m_InverseUpdated = false;
 235 }
 236
 237 // -------------------------------------------------------------------------
 238 template< class _TScalar, unsigned int _VDimension >
 239 _TScalar cpExtensions::Algorithms::
 240 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
 241 _SquaredMahalanobis( const TVector& v ) const
 242 {
 243   if( !this->m_InverseUpdated )
 244   {
 245     double d =
 246       double(
 247         vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) )
 248         );
 249     if( std::fabs( d ) > double( 0 ) )
 250       this->m_InverseUnbiasedCovariance =
 251         this->m_UnbiasedCovariance.GetInverse( );
 252     else
 253       this->m_InverseUnbiasedCovariance.SetIdentity( );
 254     this->m_InverseUpdated = true;
 255
 256   } // fi
 257   TVector x = v - this->m_Mean;
 258   return( x * ( this->m_InverseUnbiasedCovariance * x ) );
 259 }
 260
 261 #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
 262
 263 // eof - $RCSfile$