--- /dev/null
+// -------------------------------------------------------------------------
+// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
+// -------------------------------------------------------------------------
+
+#ifndef __CPEXTENSIONS__ALGORITHMS__GAUSSIANDENSITYIMAGEFILTER__H__
+#define __CPEXTENSIONS__ALGORITHMS__GAUSSIANDENSITYIMAGEFILTER__H__
+
+#include <itkImageToImageFilter.h>
+#include <itkImageScanlineConstIterator.h>
+#include <itkImageScanlineIterator.h>
+#include <itkProgressReporter.h>
+
+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$
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( );
protected:
void _AddSample( const TVector& v ) const;
_TScalar _SquaredMahalanobis( const TVector& v ) const;
+ _TScalar _Density( const TVector& v ) const;
private:
// Purposely not implemented
} // ecapseman
-#ifndef ITK_MANUAL_INSTANTIATION
#include <cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx>
-#endif // ITK_MANUAL_INSTANTIATION
#endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__
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::
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$
d #gradients=itk::CovariantVector
c cpExtensions::Algorithms::MultiScaleGaussianImageFilter< itk::Image< #scalar, #dims >, itk::Image< #gradients< #float, #dims >, #dims > >
+
+** eof
--- /dev/null
+#include <cpPluginsImageFilters/GaussianDensityImageFilter.h>
+#include <cpPlugins/Image.h>
+#include <cpExtensions/Algorithms/IterativeGaussianModelEstimator.h>
+#include <cpExtensions/Algorithms/GaussianDensityImageFilter.h>
+
+// -------------------------------------------------------------------------
+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$
--- /dev/null
+#ifndef __CPPLUGINSIMAGEFILTERS__GAUSSIANDENSITYIMAGEFILTER__H__
+#define __CPPLUGINSIMAGEFILTERS__GAUSSIANDENSITYIMAGEFILTER__H__
+
+#include <plugins/cpPluginsImageFilters/cpPluginsImageFilters_Export.h>
+#include <cpPlugins/ProcessObject.h>
+
+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$