1 // -------------------------------------------------------------------------
2 // @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
3 // -------------------------------------------------------------------------
5 #ifndef __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
6 #define __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
10 // -------------------------------------------------------------------------
11 template< class _TScalar, unsigned int _VDimension >
12 void cpExtensions::Algorithms::
13 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
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 ) );
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, ... )
33 s[ 0 ] = TScalar( x1 );
34 for( unsigned int d = 1; d < _VDimension; ++d )
35 s[ d ] = TScalar( va_arg( lst, double ) );
37 this->_AddSample( s );
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 )
48 for( unsigned d = 0; d < _VDimension; ++d )
49 s[ d ] = TScalar( array[ d ] );
50 this->_AddSample( s );
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 )
60 unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
62 for( unsigned d = 0; d < len; ++d )
63 s[ d ] = TScalar( v[ d ] );
64 this->_AddSample( s );
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 )
75 for( unsigned d = 0; d < _VDimension; ++d )
76 s[ d ] = TScalar( v[ d ] );
77 this->_AddSample( s );
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 )
88 for( unsigned d = 0; d < _VDimension; ++d )
89 s[ d ] = TScalar( p[ d ] );
90 this->_AddSample( s );
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 )
101 for( unsigned d = 0; d < _VDimension; ++d )
102 s[ d ] = TScalar( v[ d ] );
103 this->_AddSample( s );
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
116 s[ 0 ] = TScalar( x1 );
117 for( unsigned int d = 1; d < _VDimension; ++d )
118 s[ d ] = TScalar( va_arg( lst, double ) );
120 return( this->_SquaredMahalanobis( s ) );
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
131 for( unsigned d = 0; d < _VDimension; ++d )
132 s[ d ] = TScalar( array[ d ] );
133 return( this->_SquaredMahalanobis( s ) );
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
143 unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
145 for( unsigned d = 0; d < len; ++d )
146 s[ d ] = TScalar( v[ d ] );
147 return( this->_SquaredMahalanobis( s ) );
150 // -------------------------------------------------------------------------
151 template< class _TScalar, unsigned int _VDimension >
152 template< class _TOtherScalar >
153 _TScalar cpExtensions::Algorithms::
154 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
156 const itk::CovariantVector< _TOtherScalar, _VDimension >& v
160 for( unsigned d = 0; d < _VDimension; ++d )
161 s[ d ] = TScalar( v[ d ] );
162 return( this->_SquaredMahalanobis( s ) );
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
173 for( unsigned d = 0; d < _VDimension; ++d )
174 s[ d ] = TScalar( p[ d ] );
175 return( this->_SquaredMahalanobis( s ) );
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
186 for( unsigned d = 0; d < _VDimension; ++d )
187 s[ d ] = TScalar( v[ d ] );
188 return( this->_SquaredMahalanobis( s ) );
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
201 s[ 0 ] = TScalar( x1 );
202 for( unsigned int d = 1; d < _VDimension; ++d )
203 s[ d ] = TScalar( va_arg( lst, double ) );
205 return( this->_Density( s ) );
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
216 for( unsigned d = 0; d < _VDimension; ++d )
217 s[ d ] = TScalar( array[ d ] );
218 return( this->_Density( s ) );
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
228 unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension;
230 for( unsigned d = 0; d < len; ++d )
231 s[ d ] = TScalar( v[ d ] );
232 return( this->_Density( s ) );
235 // -------------------------------------------------------------------------
236 template< class _TScalar, unsigned int _VDimension >
237 template< class _TOtherScalar >
238 _TScalar cpExtensions::Algorithms::
239 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
241 const itk::CovariantVector< _TOtherScalar, _VDimension >& v
245 for( unsigned d = 0; d < _VDimension; ++d )
246 s[ d ] = TScalar( v[ d ] );
247 return( this->_Density( s ) );
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
258 for( unsigned d = 0; d < _VDimension; ++d )
259 s[ d ] = TScalar( p[ d ] );
260 return( this->_Density( s ) );
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
271 for( unsigned d = 0; d < _VDimension; ++d )
272 s[ d ] = TScalar( v[ d ] );
273 return( this->_Density( s ) );
276 // -------------------------------------------------------------------------
277 template< class _TScalar, unsigned int _VDimension >
278 cpExtensions::Algorithms::
279 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
280 IterativeGaussianModelEstimator( )
286 // -------------------------------------------------------------------------
287 template< class _TScalar, unsigned int _VDimension >
288 cpExtensions::Algorithms::
289 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
290 ~IterativeGaussianModelEstimator( )
294 // -------------------------------------------------------------------------
295 template< class _TScalar, unsigned int _VDimension >
296 void cpExtensions::Algorithms::
297 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
298 _AddSample( const TVector& v ) const
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;
305 for( unsigned int i = 0; i < _VDimension; ++i )
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 ];
312 this->m_Covariance = ( this->m_Covariance + ( m * a ) ) * b;
313 this->m_Mean = ( this->m_Mean * b ) + ( v * a );
316 TScalar( this->m_NumberOfSamples ) /
317 TScalar( this->m_NumberOfSamples - 1 );
318 this->m_UnbiasedCovariance = this->m_Covariance * u;
319 this->m_InverseUpdated = false;
322 // -------------------------------------------------------------------------
323 template< class _TScalar, unsigned int _VDimension >
324 _TScalar cpExtensions::Algorithms::
325 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
326 _SquaredMahalanobis( const TVector& v ) const
328 if( !this->m_InverseUpdated )
332 vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) )
334 if( std::fabs( d ) > double( 0 ) )
335 this->m_InverseUnbiasedCovariance =
336 this->m_UnbiasedCovariance.GetInverse( );
338 this->m_InverseUnbiasedCovariance.SetIdentity( );
339 this->m_InverseUpdated = true;
342 TVector x = v - this->m_Mean;
343 return( x * ( this->m_InverseUnbiasedCovariance * x ) );
346 // -------------------------------------------------------------------------
347 template< class _TScalar, unsigned int _VDimension >
348 _TScalar cpExtensions::Algorithms::
349 IterativeGaussianModelEstimator< _TScalar, _VDimension >::
350 _Density( const TVector& v ) const
352 static const double _2pik =
353 std::pow( double( 2 ) * double( vnl_math::pi ), _VDimension );
355 double s = double( this->_SquaredMahalanobis( v ) ) / double( 2 );
358 vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) )
360 return( _TScalar( std::exp( -s ) / std::sqrt( _2pik * d ) ) );
363 #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__