]> Creatis software - cpPlugins.git/commitdiff
...
authorLeonardo Florez-Valencia <florez-l@javeriana.edu.co>
Wed, 22 Jun 2016 03:20:30 +0000 (22:20 -0500)
committerLeonardo Florez-Valencia <florez-l@javeriana.edu.co>
Wed, 22 Jun 2016 03:20:30 +0000 (22:20 -0500)
lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h [new file with mode: 0644]
lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.h
lib/cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx
lib/cpPlugins_Instances/GaussianImageFilters.i
plugins/cpPluginsImageFilters/GaussianDensityImageFilter.cxx [new file with mode: 0644]
plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h [new file with mode: 0644]

diff --git a/lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h b/lib/cpExtensions/Algorithms/GaussianDensityImageFilter.h
new file mode 100644 (file)
index 0000000..308b3c5
--- /dev/null
@@ -0,0 +1,103 @@
+// -------------------------------------------------------------------------
+// @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$
index 95947cf0c35fe4a4ee59316c26cce522f36d8ccb..b38e06ede2d6732e86f4793a97964de10b0859b1 100644 (file)
@@ -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 <cpExtensions/Algorithms/IterativeGaussianModelEstimator.hxx>
-#endif // ITK_MANUAL_INSTANTIATION
 
 #endif // __CPEXTENSIONS__ALGORITHMS__ITERATIVEGAUSSIANMODELESTIMATOR__H__
 
index 254d201452ec6cd5ac5ad1cfa7f8215334dd4b11..acd1da76f82f380381eac997b8aed548012c4013 100644 (file)
@@ -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$
index faec2c4fe07b50675a675e72a356420937a67d77..d5b6e4188b08bd4eb2be872a613fe7ec76d1aae1 100644 (file)
@@ -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 (file)
index 0000000..12e598b
--- /dev/null
@@ -0,0 +1,80 @@
+#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$
diff --git a/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h b/plugins/cpPluginsImageFilters/GaussianDensityImageFilter.h
new file mode 100644 (file)
index 0000000..26dab40
--- /dev/null
@@ -0,0 +1,47 @@
+#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$