From 733108d258a27799df875ceda2d84c3cafea64d3 Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Tue, 21 Jun 2016 22:20:30 -0500 Subject: [PATCH] ... --- .../Algorithms/GaussianDensityImageFilter.h | 103 ++++++++++++++++++ .../IterativeGaussianModelEstimator.h | 29 ++++- .../IterativeGaussianModelEstimator.hxx | 102 +++++++++++++++++ .../GaussianImageFilters.i | 2 + .../GaussianDensityImageFilter.cxx | 80 ++++++++++++++ .../GaussianDensityImageFilter.h | 47 ++++++++ 6 files changed, 361 insertions(+), 2 deletions(-) create mode 100644 lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h create mode 100644 plugins/cpPluginsImageFilters/GaussianDensityImageFilter.cxx create mode 100644 plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h diff --git a/lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h b/lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h new file mode 100644 index 0000000..308b3c5 --- /dev/null +++ b/lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h @@ -0,0 +1,103 @@ +// ------------------------------------------------------------------------- +// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) +// ------------------------------------------------------------------------- + +#ifndef __CPEXTENSIONS__ALGORITHMS__GAUSSIANDENSITYIMAGEFILTER__H__ +#define __CPEXTENSIONS__ALGORITHMS__GAUSSIANDENSITYIMAGEFILTER__H__ + +#include +#include +#include +#include + +namespace cpExtensions +{ + namespace Algorithms + { + /** + */ + template< class _TInputImage, class _TEstimator > + class GaussianDensityImageFilter + : public itk::ImageToImageFilter< _TInputImage, itk::Image< typename _TEstimator::TScalar, _TInputImage::ImageDimension > > + { + public: + typedef typename _TEstimator::TScalar TScalar; + typedef itk::Image< TScalar, _TInputImage::ImageDimension > TOutputImage; + typedef GaussianDensityImageFilter Self; + typedef itk::ImageToImageFilter< _TInputImage, TOutputImage > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + typedef _TInputImage TInputImage; + typedef _TEstimator TEstimator; + typedef typename TOutputImage::RegionType TRegion; + + public: + itkNewMacro( Self ); + itkTypeMacro( GaussianDensityImageFilter, itkImageToImageFilter ); + + itkGetConstObjectMacro( Estimator, TEstimator ); + itkSetConstObjectMacro( Estimator, TEstimator ); + + protected: + GaussianDensityImageFilter( ) + : Superclass( ) + { + } + virtual ~GaussianDensityImageFilter( ) + { + } + + virtual void ThreadedGenerateData( + const TRegion& region, + itk::ThreadIdType threadId + ) + { + const typename TRegion::SizeType& regionSize = region.GetSize( ); + if( regionSize[ 0 ] == 0 ) + return; + const TInputImage* in = this->GetInput( ); + TOutputImage* out = this->GetOutput( 0 ); + + const size_t nLines = region.GetNumberOfPixels( ) / regionSize[ 0 ]; + itk::ProgressReporter progress( this, threadId, nLines ); + + // Define the iterators + itk::ImageScanlineConstIterator< TInputImage > inIt( in, region ); + itk::ImageScanlineIterator< TOutputImage > outIt( out, region ); + + inIt.GoToBegin( ); + outIt.GoToBegin( ); + while( !inIt.IsAtEnd( ) ) + { + while( !inIt.IsAtEndOfLine( ) ) + { + if( this->m_Estimator.IsNotNull( ) ) + outIt.Set( this->m_Estimator->Density( inIt.Get( ) ) ); + ++inIt; + ++outIt; + + } // elihw + inIt.NextLine( ); + outIt.NextLine( ); + progress.CompletedPixel( ); + + } // elihw + } + + private: + // Purposely not implemented. + GaussianDensityImageFilter( const Self& ); + void operator=( const Self& ); + + protected: + typename TEstimator::ConstPointer m_Estimator; + }; + + } // ecapseman + +} // ecapseman + +#endif // __CPEXTENSIONS__ALGORITHMS__GAUSSIANDENSITYIMAGEFILTER__H__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.h b/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.h index 95947cf..b38e06e 100644 --- a/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.h +++ b/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.h @@ -110,6 +110,32 @@ namespace cpExtensions const itk::Vector< _TOtherScalar, _VDimension >& v ) const; + template< class _TOtherScalar > + _TScalar Density( const _TOtherScalar& x1, ... ) const; + + template< class _TOtherScalar > + _TScalar Density( const _TOtherScalar* array ) const; + + template< class _TOtherScalar > + _TScalar Density( + const vnl_vector< _TOtherScalar >& v + ) const; + + template< class _TOtherScalar > + _TScalar Density( + const itk::CovariantVector< _TOtherScalar, _VDimension >& v + ) const; + + template< class _TOtherScalar > + _TScalar Density( + const itk::Point< _TOtherScalar, _VDimension >& p + ) const; + + template< class _TOtherScalar > + _TScalar Density( + const itk::Vector< _TOtherScalar, _VDimension >& v + ) const; + protected: IterativeGaussianModelEstimator( ); virtual ~IterativeGaussianModelEstimator( ); @@ -117,6 +143,7 @@ namespace cpExtensions protected: void _AddSample( const TVector& v ) const; _TScalar _SquaredMahalanobis( const TVector& v ) const; + _TScalar _Density( const TVector& v ) const; private: // Purposely not implemented @@ -137,9 +164,7 @@ namespace cpExtensions } // ecapseman -#ifndef ITK_MANUAL_INSTANTIATION #include -#endif // ITK_MANUAL_INSTANTIATION #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__ diff --git a/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx b/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx index 254d201..acd1da7 100644 --- a/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx +++ b/lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx @@ -188,6 +188,91 @@ SquaredMahalanobis( const itk::Vector< _TOtherScalar, _VDimension >& v ) const return( this->_SquaredMahalanobis( s ) ); } +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( const _TOtherScalar& x1, ... ) const +{ + TVector s; + std::va_list lst; + va_start( lst, x1 ); + s[ 0 ] = TScalar( x1 ); + for( unsigned int d = 1; d < _VDimension; ++d ) + s[ d ] = TScalar( va_arg( lst, double ) ); + va_end( lst ); + return( this->_Density( s ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( const _TOtherScalar* array ) const +{ + TVector s; + for( unsigned d = 0; d < _VDimension; ++d ) + s[ d ] = TScalar( array[ d ] ); + return( this->_Density( s ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( const vnl_vector< _TOtherScalar >& v ) const +{ + unsigned int len = ( v.size( ) < _VDimension )? v.size: _VDimension; + TVector s; + for( unsigned d = 0; d < len; ++d ) + s[ d ] = TScalar( v[ d ] ); + return( this->_Density( s ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( + const itk::CovariantVector< _TOtherScalar, _VDimension >& v + ) const +{ + TVector s; + for( unsigned d = 0; d < _VDimension; ++d ) + s[ d ] = TScalar( v[ d ] ); + return( this->_Density( s ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( const itk::Point< _TOtherScalar, _VDimension >& p ) const +{ + TVector s; + for( unsigned d = 0; d < _VDimension; ++d ) + s[ d ] = TScalar( p[ d ] ); + return( this->_Density( s ) ); +} + +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +template< class _TOtherScalar > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +Density( const itk::Vector< _TOtherScalar, _VDimension >& v ) const +{ + TVector s; + for( unsigned d = 0; d < _VDimension; ++d ) + s[ d ] = TScalar( v[ d ] ); + return( this->_Density( s ) ); +} + // ------------------------------------------------------------------------- template< class _TScalar, unsigned int _VDimension > cpExtensions::Algorithms:: @@ -258,6 +343,23 @@ _SquaredMahalanobis( const TVector& v ) const return( x * ( this->m_InverseUnbiasedCovariance * x ) ); } +// ------------------------------------------------------------------------- +template< class _TScalar, unsigned int _VDimension > +_TScalar cpExtensions::Algorithms:: +IterativeGaussianModelEstimator< _TScalar, _VDimension >:: +_Density( const TVector& v ) const +{ + static const double _2pik = + std::pow( double( 2 ) * double( vnl_math::pi ), _VDimension ); + + double s = double( this->_SquaredMahalanobis( v ) ) / double( 2 ); + double d = + double( + vnl_determinant( this->m_UnbiasedCovariance.GetVnlMatrix( ) ) + ); + return( _TScalar( std::exp( -s ) / std::sqrt( _2pik * d ) ) ); +} + #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__HXX__ // eof - $RCSfile$ diff --git a/lib/cpPlugins_Instances/GaussianImageFilters.i b/lib/cpPlugins_Instances/GaussianImageFilters.i index faec2c4..d5b6e41 100644 --- a/lib/cpPlugins_Instances/GaussianImageFilters.i +++ b/lib/cpPlugins_Instances/GaussianImageFilters.i @@ -22,3 +22,5 @@ d #scalar=#int;unsigned #int;#float d #gradients=itk::CovariantVector c cpExtensions::Algorithms::MultiScaleGaussianImageFilter< itk::Image< #scalar, #dims >, itk::Image< #gradients< #float, #dims >, #dims > > + +** eof diff --git a/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.cxx b/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.cxx new file mode 100644 index 0000000..12e598b --- /dev/null +++ b/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.cxx @@ -0,0 +1,80 @@ +#include +#include +#include +#include + +// ------------------------------------------------------------------------- +cpPluginsImageFilters::GaussianDensityImageFilter:: +GaussianDensityImageFilter( ) + : Superclass( ) +{ + this->_AddInput( "Image" ); + this->_AddInput( "Estimator" ); + this->_AddOutput< cpPlugins::Image >( "Output" ); + + /* TODO + this->m_Parameters.ConfigureAsReal( "LowerThresholdValue" ); + this->m_Parameters.ConfigureAsReal( "UpperThresholdValue" ); + this->m_Parameters.ConfigureAsUint( "InsideValue" ); + this->m_Parameters.ConfigureAsUint( "OutsideValue" ); + + this->m_Parameters.SetReal( "LowerThresholdValue", 0 ); + this->m_Parameters.SetReal( "UpperThresholdValue", 10000 ); + this->m_Parameters.SetReal( "InsideValue", 1 ); + this->m_Parameters.SetReal( "OutsideValue", 0 ); + this->m_Parameters.SetSelectedChoice( "OutputResolution", "unsigned char" ); + */ +} + +// ------------------------------------------------------------------------- +cpPluginsImageFilters::GaussianDensityImageFilter:: +~GaussianDensityImageFilter( ) +{ +} + +// ------------------------------------------------------------------------- +void cpPluginsImageFilters::GaussianDensityImageFilter:: +_GenerateData( ) +{ + auto image = this->GetInputData< itk::DataObject >( "Image" ); + cpPlugins_Image_Demangle_Pixel_AllScalars ( _GD0, image, 2 ); + else cpPlugins_Image_Demangle_Pixel_AllScalars( _GD0, image, 3 ); + else this->_Error( "No valid input image." ); +} + +// ------------------------------------------------------------------------- +template< class _TImage > +void cpPluginsImageFilters::GaussianDensityImageFilter:: +_GD0( _TImage* image ) +{ + static const unsigned int Dim = 1; + typedef cpExtensions::Algorithms::IterativeGaussianModelEstimator< float, Dim > _TFloatEst; + typedef cpExtensions::Algorithms::IterativeGaussianModelEstimator< double, Dim > _TDoubleEst; + + auto float_est = this->GetInputData< _TFloatEst >( "Estimator" ); + auto double_est = this->GetInputData< _TDoubleEst >( "Estimator" ); + if( float_est != NULL ) + this->_GD1( image, float_est ); + else if( double_est != NULL ) + this->_GD1( image, double_est ); + else + this->_Error( "Invalid gaussian model estimator." ); +} + +// ------------------------------------------------------------------------- +template< class _TImage, class _TEstimator > +void cpPluginsImageFilters::GaussianDensityImageFilter:: +_GD1( _TImage* image, _TEstimator* estimator ) +{ + typedef + cpExtensions::Algorithms::GaussianDensityImageFilter< _TImage, _TEstimator > + _TFilter; + + auto filter = this->_CreateITK< _TFilter >( ); + filter->SetInput( image ); + filter->SetEstimator( estimator ); + filter->Update( ); + this->GetOutput( "Output" )->SetITK( filter->GetOutput( ) ); +} + +// eof - $RCSfile$ diff --git a/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h b/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h new file mode 100644 index 0000000..26dab40 --- /dev/null +++ b/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h @@ -0,0 +1,47 @@ +#ifndef __CPPLUGINSIMAGEFILTERS__GAUSSIANDENSITYIMAGEFILTER__H__ +#define __CPPLUGINSIMAGEFILTERS__GAUSSIANDENSITYIMAGEFILTER__H__ + +#include +#include + +namespace cpPluginsImageFilters +{ + /** + */ + class cpPluginsImageFilters_EXPORT GaussianDensityImageFilter + : public cpPlugins::ProcessObject + { + public: + typedef GaussianDensityImageFilter Self; + typedef cpPlugins::ProcessObject Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + public: + itkNewMacro( Self ); + itkTypeMacro( GaussianDensityImageFilter, cpPlugins::ProcessObject ); + cpPlugins_Id_Macro( GaussianDensityImageFilter, ImageFilters ); + + protected: + GaussianDensityImageFilter( ); + virtual ~GaussianDensityImageFilter( ); + + virtual void _GenerateData( ) ITK_OVERRIDE; + + template< class _TImage > + inline void _GD0( _TImage* image ); + + template< class _TImage, class _TEstimator > + inline void _GD1( _TImage* image, _TEstimator* estimator ); + + private: + // Purposely not implemented + GaussianDensityImageFilter( const Self& ); + Self& operator=( const Self& ); + }; + +} // ecapseman + +#endif // __CPPLUGINSIMAGEFILTERS__GAUSSIANDENSITYIMAGEFILTER__H__ + +// eof - $RCSfile$ -- 2.45.1