From: Leonardo Florez-Valencia Date: Tue, 28 Apr 2015 22:30:30 +0000 (-0500) Subject: Bezier function added. X-Git-Tag: v0.1~401 X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=commitdiff_plain;h=263279e49b0b075a6944ddd1c498d87ef817e118;p=cpPlugins.git Bezier function added. --- diff --git a/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.h b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.h new file mode 100644 index 0000000..641220c --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.h @@ -0,0 +1,78 @@ +// ------------------------------------------------------------------------- +// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +// ------------------------------------------------------------------------- + +#ifndef __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__ +#define __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__ + +#include +#include +#include + +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 + +#endif // __CPPLUGINS__EXTENSIONS__ALGORITHMS__BEZIERCURVEFUNCTION__H__ + +// eof - $RCSfile$ diff --git a/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx new file mode 100644 index 0000000..dad1482 --- /dev/null +++ b/lib/cpPlugins/Extensions/Algorithms/BezierCurveFunction.hxx @@ -0,0 +1,194 @@ +// ------------------------------------------------------------------------- +// @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$ diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h index 94e5495..017691e 100644 --- a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.h @@ -47,14 +47,20 @@ namespace cpPlugins 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; @@ -75,24 +81,23 @@ namespace cpPlugins 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 diff --git a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx index bf87054..a8006a0 100644 --- a/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx +++ b/lib/cpPlugins/Extensions/Algorithms/IterativeGaussianModelEstimator.hxx @@ -7,6 +7,7 @@ #include #include +#include #include #include @@ -14,65 +15,95 @@ 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 ); } // ------------------------------------------------------------------------- @@ -81,13 +112,16 @@ template< class V > 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 ) ) ); @@ -95,9 +129,10 @@ Probability( const V& sample ) const 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 @@ -127,9 +162,11 @@ void 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 ] ); @@ -142,12 +179,12 @@ void 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( ); } @@ -162,9 +199,24 @@ AddSample( const V& sample ) 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( ); @@ -204,6 +256,42 @@ cpPlugins::Extensions::Algorithms::IterativeGaussianModelEstimator< S, D >:: { } +// ------------------------------------------------------------------------- +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$ diff --git a/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h b/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h index 9563a63..56ffd5c 100644 --- a/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h +++ b/lib/cpPlugins/Extensions/Algorithms/RGBToYPbPrFunction.h @@ -33,9 +33,9 @@ namespace cpPlugins { 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 );