--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__
+#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__
+
+#include <vector>
+#include <itkFunctionBase.h>
+#include <itkMatrix.h>
+
+namespace cpPlugins
+{
+ namespace Extensions
+ {
+ namespace Algorithms
+ {
+ /**
+ * Uses the De Casteljau's algorithm.
+ */
+ template< class V >
+ class BezierCurveFunction
+ : public itk::FunctionBase< typename V::ValueType, V >
+ {
+ public:
+ typedef BezierCurveFunction Self;
+ typedef itk::FunctionBase< typename V::ValueType, V > Superclass;
+ typedef itk::SmartPointer< Self > Pointer;
+ typedef itk::SmartPointer< const Self > ConstPointer;
+
+ typedef V TVector;
+ typedef typename V::ValueType TScalar;
+
+ typedef
+ itk::Matrix< TScalar, TVector::Dimension, TVector::Dimension >
+ TFrame;
+ typedef std::vector< TVector > TVectorsContainer;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( BezierCurveFunction, itkFunctionBase );
+
+ public:
+ virtual void AddPoint( const TVector& v );
+ virtual unsigned int GetNumberOfPoints( ) const;
+
+ virtual TVector Evaluate( const TScalar& u ) const;
+ virtual TFrame EvaluateFrenetFrame( const TScalar& u ) const;
+ virtual TScalar EvaluateLength( ) const;
+
+ protected:
+ BezierCurveFunction( );
+ virtual ~BezierCurveFunction( ) { }
+
+ void _UpdateDerivative( ) const;
+
+ private:
+ // Purposely not implemented
+ BezierCurveFunction( const Self& other );
+ Self& operator=( const Self& other );
+
+ protected:
+ TVectorsContainer m_Vectors;
+ mutable Pointer m_Derivative;
+ mutable bool m_DerivativeUpdated;
+ };
+
+ } // ecapseman
+
+ } // ecapseman
+
+} // ecapseman
+
+#include <cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx>
+
+#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__
+
+// eof - $RCSfile$
--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
+#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
+
+// -------------------------------------------------------------------------
+template< class V >
+void cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+AddPoint( const TVector& v )
+{
+ this->m_Vectors.push_back( v );
+ this->m_DerivativeUpdated = false;
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+unsigned int cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+GetNumberOfPoints( ) const
+{
+ return( this->m_Vectors.size( ) );
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+TVector cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+Evaluate( const TScalar& u ) const
+{
+ TVectorsContainer Q = this->m_Vectors;
+ unsigned int n = Q.size( );
+ TScalar _1u = TScalar( 1 ) - u;
+
+ for( unsigned int k = 1; k < n; k++ )
+ {
+ // CM Fixed a bug appearing under Windows : changed the stopping condition from <= to <.
+ // Otherwise, on the last step, an element out of the range of vector Q is accessed (Q[ i + 1 ])...
+ for( unsigned int i = 0; i < n - k; i++ )
+ Q[ i ] = ( Q[ i ] * _1u ) + ( Q[ i + 1 ] * u );
+
+ } // rof
+ return( Q[ 0 ] );
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+TFrame cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+EvaluateFrenetFrame( const TScalar& u ) const
+{
+ TFrame fr;
+ fr.Fill( TScalar( 0 ) );
+ if( TVector::Dimension == 2 )
+ {
+ this->_UpdateDerivative( );
+ this->m_Derivative->_UpdateDerivative( );
+
+ TVector vT = this->m_Derivative->Evaluate( u );
+ TScalar nvT = vT.GetNorm( );
+ if( TScalar( 0 ) < nvT )
+ vT /= nvT;
+
+ fr[ 0 ][ 0 ] = vT[ 0 ];
+ fr[ 1 ][ 0 ] = vT[ 1 ];
+
+ fr[ 0 ][ 1 ] = -vT[ 1 ];
+ fr[ 1 ][ 1 ] = vT[ 0 ];
+
+ } // fi
+
+ /* TODO
+ if( TVector::Dimension == 3 )
+ {
+ this->_UpdateDerivative( );
+ this->m_Derivative->_UpdateDerivative( );
+
+ TVector vT = this->m_Derivative->Evaluate( u );
+ TVector vN = this->m_Derivative->m_Derivative->Evaluate( u );
+ TScalar nvT = vT.GetNorm( );
+ TScalar nvN = vN.GetNorm( );
+
+ if( nvT > TScalar( 0 ) && nvN > TScalar( 0 ) )
+ {
+ vT /= nvT;
+ vN /= nvN;
+
+ TVector vB;
+ vB[ 0 ] = ( vT[ 1 ] * vN[ 2 ] ) - ( vT[ 2 ] * vN[ 1 ] );
+ vB[ 1 ] = ( vT[ 2 ] * vN[ 0 ] ) - ( vT[ 0 ] * vN[ 2 ] );
+ vB[ 2 ] = ( vT[ 0 ] * vN[ 1 ] ) - ( vT[ 1 ] * vN[ 0 ] );
+
+ TScalar nvB = vB.GetNorm( );
+ if( nvB > TScalar( 0 ) )
+ {
+ vB /= nvB;
+
+ // WARNING: Paranoiac test
+ vN[ 0 ] = ( vB[ 1 ] * vT[ 2 ] ) - ( vB[ 2 ] * vT[ 1 ] );
+ vN[ 1 ] = ( vB[ 2 ] * vT[ 0 ] ) - ( vB[ 0 ] * vT[ 2 ] );
+ vN[ 2 ] = ( vB[ 0 ] * vT[ 1 ] ) - ( vB[ 1 ] * vT[ 0 ] );
+
+ for( unsigned int d = 0; d < 3; d++ )
+ {
+ fr[ d ][ 0 ] = vT[ d ];
+ fr[ d ][ 1 ] = vN[ d ];
+ fr[ d ][ 2 ] = vB[ d ];
+
+ } // rof
+ }
+ else
+ {
+ // WARNING: Trick to avoid numerical instabilities
+ // in straight lines.
+ typedef itk::Vector< TScalar, 3 > _TVector3;
+ typedef ext::VectorToFrameFunction< _TVector3 > _TFunction;
+
+ _TVector3 vT3;
+ vT3[ 0 ] = vT[ 0 ];
+ vT3[ 1 ] = vT[ 1 ];
+ vT3[ 2 ] = vT[ 2 ];
+
+ typename _TFunction::Pointer fun = _TFunction::New( );
+ typename _TFunction::TFrame ffr = fun->Evaluate( vT3 );
+
+ fr[ 0 ][ 0 ] = ffr[ 0 ][ 0 ];
+ fr[ 0 ][ 1 ] = ffr[ 0 ][ 1 ];
+ fr[ 0 ][ 2 ] = ffr[ 0 ][ 2 ];
+
+ fr[ 1 ][ 0 ] = ffr[ 1 ][ 0 ];
+ fr[ 1 ][ 1 ] = ffr[ 1 ][ 1 ];
+ fr[ 1 ][ 2 ] = ffr[ 1 ][ 2 ];
+
+ fr[ 2 ][ 0 ] = ffr[ 2 ][ 0 ];
+ fr[ 2 ][ 1 ] = ffr[ 2 ][ 1 ];
+ fr[ 2 ][ 2 ] = ffr[ 2 ][ 2 ];
+
+ } // fi
+
+ } // fi
+
+ } // fi
+ */
+ return( fr );
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+typename cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+TScalar cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+EvaluateLength( ) const
+{
+ unsigned int n = this->GetNumberOfPoints( ) << 1;
+ TScalar d = TScalar( 0 );
+ TVector v0 = this->Evaluate( 0 );
+ for( unsigned int i = 1; i < n; i++ )
+ {
+ TVector v1 = this->Evaluate( TScalar( i ) / TScalar( n - 1 ) );
+ d += ( v1 - v0 ).GetNorm( );
+ v0 = v1;
+
+ } // rof
+ return( d );
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+BezierCurveFunction( )
+ : Superclass( ),
+ m_DerivativeUpdated( false )
+{
+}
+
+// -------------------------------------------------------------------------
+template< class V >
+void cpPlugins::Extensions::Algorithms::BezierCurveFunction< V >::
+_UpdateDerivative( ) const
+{
+ if( this->m_DerivativeUpdated )
+ return;
+
+ this->m_Derivative = Self::New( );
+ unsigned int n = this->m_Vectors.size( ) - 1;
+ for( unsigned int i = 0; i < n; i++ )
+ this->m_Derivative->AddPoint(
+ TScalar( n ) * ( this->m_Vectors[ i + 1 ] - this->m_Vectors[ i ] )
+ );
+}
+
+#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__HXX__
+
+// eof - $RCSfile$
itkNewMacro( Self );
itkTypeMacro( IterativeGaussianModelEstimator, itkObject );
- itkGetConstMacro( MinimumProbability, S );
- itkGetConstMacro( MaximumProbability, S );
-
public:
unsigned long GetNumberOfSamples( ) const
- { return( this->m_Samples.size( ) ); }
+ { return( ( unsigned long )( this->m_N ) ); }
+ const TMatrix& GetMu( ) const
+ { return( this->m_M ); }
+ const TMatrix& GetOmega( ) const
+ { return( this->m_O ); }
+
+ void SetNumberOfSamples( unsigned long n );
+ void SetMu( const TMatrix& m );
+ void SetOmega( const TMatrix& O );
- void UpdateModel( );
+ bool SaveModelToFile( const std::string& filename ) const;
+ bool LoadModelFromFile( const std::string& filename );
template< class V >
S Probability( const V& sample ) const;
IterativeGaussianModelEstimator( );
virtual ~IterativeGaussianModelEstimator( );
+ protected:
+ void _UpdateModel( ) const;
+
private:
// Purposely not implemented
IterativeGaussianModelEstimator( const Self& other );
void operator=( const Self& other );
protected:
- TMatrix m_S1;
- TMatrix m_S2;
- TMatrices m_Samples;
-
- bool m_Updated;
- TMatrix m_Cov;
- TMatrix m_Mean;
- TMatrix m_InvCov;
- S m_DensityCoeff;
-
- S m_MinimumProbability;
- S m_MaximumProbability;
+ S m_N;
+ TMatrix m_M;
+ TMatrix m_O;
+
+ mutable bool m_Updated;
+ mutable TMatrix m_Cov;
+ mutable TMatrix m_Inv;
+ mutable S m_Norm;
};
} // ecapseman
#include <cmath>
#include <cstdarg>
+#include <fstream>
#include <vnl/vnl_math.h>
#include <vnl/vnl_inverse.h>
template< class S, unsigned int D >
void
cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
-UpdateModel( )
+SetNumberOfSamples( unsigned long n )
{
- if( this->m_Updated )
- return;
-
- this->m_Cov.set_size( D, D );
- this->m_Mean.set_size( D, 1 );
-
- // Compute covariance matrix and mean vector
- unsigned long N = this->m_Samples.size( );
- this->m_Mean = this->m_S1;
- if( N > 0 )
- this->m_Mean /= S( N );
- if( N > 1 )
+ this->m_N = S( n );
+ this->m_Updated = false;
+ this->Modified( );
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void
+cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
+SetMu( const TMatrix& m )
+{
+ if( m.rows( ) == D && m.columns( ) == 1 )
{
- this->m_Cov = this->m_S2;
- this->m_Cov -= ( this->m_S1 * this->m_S1.transpose( ) ) / S( N );
- this->m_Cov /= S( N - 1 );
+ this->m_M = m;
+ this->m_Updated = false;
+ this->Modified( );
}
else
- this->m_Cov.fill( S( 0 ) );
+ {
+ itkExceptionMacro(
+ << "Input Mu matrix is not a " << D << "x1 matrix"
+ );
- // Compute inverse and determinant
- S det = vnl_determinant( this->m_Cov );
- if( !( det > S( 0 ) ) )
+ } // fi
+}
+
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void
+cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
+SetOmega( const TMatrix& O )
+{
+ if( O.rows( ) == D && O.columns( ) == D )
{
- this->m_InvCov.set_size( D, D );
- this->m_InvCov.fill( S( 0 ) );
- this->m_DensityCoeff = S( 0 );
+ this->m_O = O;
+ this->m_Updated = false;
+ this->Modified( );
}
else
{
- this->m_InvCov = vnl_inverse( this->m_Cov );
- double _2piD = std::pow( double( 2 ) * double( vnl_math::pi ), D );
- this->m_DensityCoeff = S( 1 ) / S( std::sqrt( _2piD * double( det ) ) );
+ itkExceptionMacro(
+ << "Input Omega matrix is not a " << D << "x" << D << " matrix"
+ );
} // fi
+}
- // Compute minimum and maximum probabilities from input samples
- static S sample[ D ];
- for( unsigned long i = 0; i < this->m_Samples.size( ); ++i )
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+bool
+cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
+SaveModelToFile( const std::string& filename ) const
+{
+ std::ofstream out( filename.c_str( ) );
+ if( out )
{
- for( unsigned int d = 0; d < D; ++d )
- sample[ d ] = this->m_Samples[ i ][ d ][ 0 ];
- S p = this->Probability( sample );
- if( i == 0 )
- {
- this->m_MinimumProbability = p;
- this->m_MaximumProbability = p;
- }
- else
- {
- if( p < this->m_MinimumProbability ) this->m_MinimumProbability = p;
- if( this->m_MaximumProbability < p ) this->m_MaximumProbability = p;
-
- } // fi
+ out << this->m_N << std::endl;
+ out << this->m_M << std::endl;
+ out << this->m_O << std::endl;
+ out.close( );
+ return( true );
+ }
+ else
+ return( false );
+}
- } // rof
- this->m_Updated = true;
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+bool
+cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
+LoadModelFromFile( const std::string& filename )
+{
+ std::ifstream in( filename.c_str( ) );
+ if( in )
+ {
+ this->Clear( );
+ in >> this->m_N;
+ for( unsigned int i = 0; i < D; ++i )
+ in >> this->m_M[ i ][ 0 ];
+ for( unsigned int j = 0; j < D; ++j )
+ for( unsigned int i = 0; i < D; ++i )
+ in >> this->m_O[ i ][ j ];
+ in.close( );
+ return( true );
+ }
+ else
+ return( false );
}
// -------------------------------------------------------------------------
S cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
Probability( const V& sample ) const
{
+ if( !this->m_Updated )
+ this->_UpdateModel( );
+
TMatrix c( D, 1 );
for( unsigned int d = 0; d < D; ++d )
- c[ d ][ 0 ] = S( sample[ d ] ) - this->m_Mean[ d ][ 0 ];
- if( S( 0 ) < this->m_DensityCoeff )
+ c[ d ][ 0 ] = S( sample[ d ] ) - this->m_M[ d ][ 0 ];
+ if( S( 0 ) < this->m_Norm )
{
// Covariance is NOT null
- double v = double( ( c.transpose( ) * ( this->m_InvCov * c ) )[ 0 ][ 0 ] );
+ double v = double( ( c.transpose( ) * ( this->m_Inv * c ) )[ 0 ][ 0 ] );
v /= double( 2 );
return( S( std::exp( -v ) ) );
else
{
// Covariance is null
- bool equal = true;
- for( unsigned int d = 0; d < D && equal; ++d )
- equal &= !( S( 0 ) < S( std::fabs( double( c[ d ][ 0 ] ) ) ) );
+ S n = S( 0 );
+ for( unsigned int d = 0; d < D; ++d )
+ n += c[ d ][ 0 ] * c[ d ][ 0 ];
+ bool equal = ( double( n ) < double( 1e-10 ) );
return( ( equal )? S( 1 ): S( 0 ) );
} // fi
cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
GetModel( V& m, M& E ) const
{
+ if( !this->m_Updated )
+ this->_UpdateModel( );
for( unsigned int i = 0; i < D; ++i )
{
- m[ i ] = double( this->m_Mean[ i ][ 0 ] );
+ m[ i ] = double( this->m_M[ i ][ 0 ] );
for( unsigned int j = 0; j < D; ++j )
E[ i ][ j ] = double( this->m_Cov[ i ][ j ] );
cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
Clear( )
{
+ this->m_N = S( 0 );
+ this->m_M.set_size( D, 1 );
+ this->m_O.set_size( D, D );
+ this->m_M.fill( S( 0 ) );
+ this->m_O.fill( S( 0 ) );
this->m_Updated = false;
- this->m_Samples.clear( );
- this->m_S1.set_size( D, 1 );
- this->m_S2.set_size( D, D );
- this->m_S1.fill( S( 0 ) );
- this->m_S2.fill( S( 0 ) );
this->Modified( );
}
for( unsigned int d = 0; d < D; ++d )
s[ d ][ 0 ] = S( sample[ d ] );
- this->m_Samples.push_back( s );
- this->m_S1 += s;
- this->m_S2 += s * s.transpose( );
+ // Update mean
+ S coeff = this->m_N;
+ this->m_N += S( 1 );
+ coeff /= this->m_N;
+ this->m_M = ( this->m_M * coeff ) + ( s / this->m_N );
+
+ // Update omega operand
+ if( this->m_N == 1 )
+ this->m_O = s * s.transpose( );
+ else if( this->m_N == 2 )
+ this->m_O += s * s.transpose( );
+ else
+ {
+ S inv = S( 1 ) / ( this->m_N - S( 1 ) );
+ this->m_O = this->m_O * ( this->m_N - S( 2 ) ) * inv;
+ this->m_O += ( s * s.transpose( ) ) * inv;
+
+ } // fi
this->m_Updated = false;
this->Modified( );
{
}
+// -------------------------------------------------------------------------
+template< class S, unsigned int D >
+void
+cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >::
+_UpdateModel( ) const
+{
+ static const double _2piD =
+ std::pow( double( 2 ) * double( vnl_math::pi ), D );
+
+ // Update covariance
+ this->m_Cov =
+ this->m_O -
+ (
+ ( this->m_M * this->m_M.transpose( ) ) *
+ ( this->m_N / ( this->m_N - S( 1 ) ) )
+ );
+
+ // Compute inverse and determinant
+ S det = vnl_determinant( this->m_Cov );
+ if( !( det > S( 0 ) ) )
+ {
+ this->m_Inv.set_size( D, D );
+ this->m_Inv.fill( S( 0 ) );
+ this->m_Norm = S( 0 );
+ }
+ else
+ {
+ this->m_Inv = vnl_inverse( this->m_Cov );
+ this->m_Norm = S( 1 ) / S( std::sqrt( _2piD * double( det ) ) );
+
+ } // fi
+
+ // Object now is updated
+ this->m_Updated = true;
+}
+
#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__
// eof - $RCSfile$
{
static const O M[] =
{
- O( 0.299 ), O( 0.587 ), O( 0.114 ),
- O( -0.169 ), O( -0.331 ), O( 0.500 ),
- O( 0.500 ), O( -0.419 ), O( -0.081 )
+ O( 0.2990 ), O( 0.5870 ), O( 0.1140 ),
+ O( -0.1687 ), O( -0.3313 ), O( 0.5000 ),
+ O( 0.5000 ), O( -0.4187 ), O( -0.0813 )
};
static const vnl_matrix< O > vM( M, 3, 3 );
static const itk::Matrix< O, 3, 3 > iM( vM );
git://git.creatis.insa-lyon.fr / cpPlugins.git / commitdiff