From 1df9496bf2f5d456616947d377f991a1022ebf6f Mon Sep 17 00:00:00 2001 From: Leonardo Florez-Valencia Date: Tue, 21 Jun 2016 17:27:27 -0500 Subject: [PATCH] Even more plugins. --- lib/cpExtensions/Algorithms/FluxMedialness.h | 74 ++++++ .../Algorithms/FluxMedialness.hxx | 108 ++++++++ .../Algorithms/GradientFunctionBase.h | 84 ------- .../Algorithms/GradientFunctionBase.hxx | 82 ------ .../Algorithms/GradientImageFunctionBase.h | 68 +++++ .../Algorithms/GradientImageFunctionBase.hxx | 59 +++++ .../Algorithms/GulsunTekMedialness.h | 53 ++-- .../Algorithms/GulsunTekMedialness.hxx | 235 ++++++------------ lib/cpExtensions/Algorithms/MFluxMedialness.h | 74 ++++++ .../Algorithms/MFluxMedialness.hxx | 159 ++++++++++++ 10 files changed, 640 insertions(+), 356 deletions(-) create mode 100644 lib/cpExtensions/Algorithms/FluxMedialness.h create mode 100644 lib/cpExtensions/Algorithms/FluxMedialness.hxx delete mode 100644 lib/cpExtensions/Algorithms/GradientFunctionBase.h delete mode 100644 lib/cpExtensions/Algorithms/GradientFunctionBase.hxx create mode 100644 lib/cpExtensions/Algorithms/GradientImageFunctionBase.h create mode 100644 lib/cpExtensions/Algorithms/GradientImageFunctionBase.hxx create mode 100644 lib/cpExtensions/Algorithms/MFluxMedialness.h create mode 100644 lib/cpExtensions/Algorithms/MFluxMedialness.hxx diff --git a/lib/cpExtensions/Algorithms/FluxMedialness.h b/lib/cpExtensions/Algorithms/FluxMedialness.h new file mode 100644 index 0000000..f71da00 --- /dev/null +++ b/lib/cpExtensions/Algorithms/FluxMedialness.h @@ -0,0 +1,74 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__ +#define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__ + +#include + +namespace cpExtensions +{ + namespace Algorithms + { + /** + */ + template< class _TGradient > + class FluxMedialness + : public GradientImageFunctionBase< _TGradient > + { + public: + typedef FluxMedialness Self; + typedef GradientImageFunctionBase< _TGradient > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + itkStaticConstMacro( Dimension, unsigned int, Superclass::Dimension ); + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TScalar TScalar; + typedef typename Superclass::TIndex TIndex; + typedef typename Superclass::TVector TVector; + typedef typename Superclass::TPoint TPoint; + + typedef std::vector< double > TRCandidates; + + public: + itkNewMacro( Self ); + itkTypeMacro( FluxMedialness, GradientImageFunctionBase ); + + itkGetConstMacro( RadiusStep, double ); + itkGetConstMacro( MinRadius, double ); + itkGetConstMacro( MaxRadius, double ); + itkGetConstMacro( RadialSampling, unsigned int ); + + itkSetMacro( RadiusStep, double ); + itkSetMacro( MinRadius, double ); + itkSetMacro( MaxRadius, double ); + itkSetMacro( RadialSampling, unsigned int ); + + protected: + FluxMedialness( ); + virtual ~FluxMedialness( ); + + virtual TOutput _Evaluate( const TIndex& i ) const ITK_OVERRIDE; + + private: + // Purposely not implemented. + FluxMedialness( const Self& ); + void operator=( const Self& ); + + protected: + double m_MinRadius; + double m_MaxRadius; + unsigned int m_RadialSampling; + double m_RadiusStep; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/FluxMedialness.hxx b/lib/cpExtensions/Algorithms/FluxMedialness.hxx new file mode 100644 index 0000000..f4d8b9c --- /dev/null +++ b/lib/cpExtensions/Algorithms/FluxMedialness.hxx @@ -0,0 +1,108 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__ +#define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::FluxMedialness< _TGradient >:: +FluxMedialness( ) + : Superclass( ), + m_MinRadius( double( 0 ) ), + m_MaxRadius( double( 1 ) ), + m_RadialSampling( 4 ), + m_RadiusStep( double( 1 ) ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::FluxMedialness< _TGradient >:: +~FluxMedialness( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +typename cpExtensions::Algorithms::FluxMedialness< _TGradient >:: +TOutput cpExtensions::Algorithms::FluxMedialness< _TGradient >:: +_Evaluate( const TIndex& i ) const +{ + itk::Object::GlobalWarningDisplayOff( ); + + double pi2n = double( 2 ) * double( vnl_math::pi ); + pi2n /= int( this->m_RadialSampling ); + const _TGradient* img = this->GetInputImage( ); + // Gradient in central pixel + //TVector grad_idx = img->GetPixel( i ); + + double Flux1 = 0; + double Flux = 0; + + TRCandidates FluxFinal; + TRCandidates radiusGenerated; + double dR = double( 0 ); + double optR = double( 0 ); + TPoint center; + img->TransformIndexToPhysicalPoint( i, center ); + double radius = double(0); + + for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ ) + { + for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ ) + { + dR = double( 0 ); + FluxFinal.clear(); + radiusGenerated.clear(); + radius = this->m_MinRadius; + while( radius <= this->m_MaxRadius ) + { + Flux = 0; + for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling ; I_radial++ ) + { + Flux1 = 0; + + // Direction of first profile + typename TPoint::VectorType dir1; + dir1.Fill( double( 0 ) ); + dir1[ cx ] = std::cos( pi2n * double( I_radial ) ); + dir1[ cy ] = std::sin( pi2n * double( I_radial ) ); + dir1 *= (radius); + TIndex rIdx; + if (img->TransformPhysicalPointToIndex( center + dir1, rIdx )) + { + TVector grad_rIdx = img->GetPixel( rIdx ); + TVector u_i1; + u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) ); + u_i1.Normalize( ); + // dot product + Flux1 = grad_rIdx * u_i1; + } + + Flux += Flux1; + + } // rof + //std::cout<<" radius:"<m_RadialSampling; + FluxFinal.push_back(Flux); + radiusGenerated.push_back(radius); + radius += this->m_RadiusStep; + + } //elihw + + dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) ); + optR= (dR>optR)? dR:optR; + + } // rof + + } // rof + return( TScalar(optR) ); +} + +#endif // __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/GradientFunctionBase.h b/lib/cpExtensions/Algorithms/GradientFunctionBase.h deleted file mode 100644 index 0472cc4..0000000 --- a/lib/cpExtensions/Algorithms/GradientFunctionBase.h +++ /dev/null @@ -1,84 +0,0 @@ -// ------------------------------------------------------------------------- -// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) -// ------------------------------------------------------------------------- - -#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__ -#define __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__ - -#include -#include - -namespace cpExtensions -{ - namespace Algorithms - { - /** - */ - template< class G > - class GradientFunctionBase - : public itk::ImageFunction< G, typename G::PixelType::ValueType, typename G::PixelType::ValueType > - { - public: - // Types from input arguments - typedef G TGradient; - typedef typename G::PixelType TVector; - typedef typename TVector::ValueType TScalar; - itkStaticConstMacro( Dimension, unsigned int, G::ImageDimension ); - - // Standard itk types - typedef GradientFunctionBase Self; - typedef itk::ImageFunction< G, TScalar, TScalar > Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - // Types from base itk::ImageFunction - typedef typename Superclass::InputType TInput; - typedef typename Superclass::OutputType TOutput; - typedef typename Superclass::PointType TPoint; - typedef typename Superclass::ContinuousIndexType TContIndex; - typedef typename Superclass::IndexType TIndex; - - // Sparse buffer - typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer; - - public: - itkTypeMacro( GradientFunctionBase, itkImageFunction ); - - itkBooleanMacro( BufferResults ); - itkGetConstMacro( BufferResults, bool ); - itkSetMacro( BufferResults, bool ); - - public: - virtual void ResetBuffer( ); - - virtual TOutput Evaluate( const TPoint& p ) const; - virtual TOutput EvaluateAtIndex( const TIndex& i ) const; - virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const; - - protected: - GradientFunctionBase( ); - virtual ~GradientFunctionBase( ); - - virtual TOutput _Evaluate( const TIndex& i ) const = 0; - - private: - // Purposely not implemented. - GradientFunctionBase( const Self& ); - void operator=( const Self& ); - - protected: - mutable TBuffer m_Buffer; - bool m_BufferResults; - }; - - } // ecapseman - -} // ecapseman - -#ifndef ITK_MANUAL_INSTANTIATION -#include -#endif // ITK_MANUAL_INSTANTIATION - -#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__ - -// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/GradientFunctionBase.hxx b/lib/cpExtensions/Algorithms/GradientFunctionBase.hxx deleted file mode 100644 index cc0bac3..0000000 --- a/lib/cpExtensions/Algorithms/GradientFunctionBase.hxx +++ /dev/null @@ -1,82 +0,0 @@ -// ------------------------------------------------------------------------- -// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) -// ------------------------------------------------------------------------- - -#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__HXX__ -#define __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__HXX__ - -// ------------------------------------------------------------------------- -template< class G > -void cpExtensions::Algorithms::GradientFunctionBase< G >:: -ResetBuffer( ) -{ - this->m_Buffer.clear( ); -} - -// ------------------------------------------------------------------------- -template< class G > -typename cpExtensions::Algorithms::GradientFunctionBase< G >:: -TOutput cpExtensions::Algorithms::GradientFunctionBase< G >:: -Evaluate( const TPoint& p ) const -{ - TIndex i; - this->GetInputImage( )->TransformPhysicalPointToIndex( p, i ); - return( this->EvaluateAtIndex( i ) ); -} - -// ------------------------------------------------------------------------- -template< class G > -typename cpExtensions::Algorithms::GradientFunctionBase< G >:: -TOutput cpExtensions::Algorithms::GradientFunctionBase< G >:: -EvaluateAtIndex( const TIndex& i ) const -{ - TOutput res = TOutput( 0 ); - bool computed = false; - if( this->m_BufferResults ) - { - typename TBuffer::const_iterator bIt = this->m_Buffer.find( i ); - computed = ( bIt != this->m_Buffer.end( ) ); - res = ( computed )? bIt->second: res; - - } // fi - - if( !computed ) - res = this->_Evaluate( i ); - - if( this->m_BufferResults ) - this->m_Buffer[ i ] = res; - return( res ); -} - -// ------------------------------------------------------------------------- -template< class G > -typename cpExtensions::Algorithms::GradientFunctionBase< G >:: -TOutput cpExtensions::Algorithms::GradientFunctionBase< G >:: -EvaluateAtContinuousIndex( const TContIndex& i ) const -{ - TPoint p; - this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p ); - return( this->Evaluate( p ) ); -} - -// ------------------------------------------------------------------------- -template< class G > -cpExtensions::Algorithms::GradientFunctionBase< G >:: -GradientFunctionBase( ) - : Superclass( ), - m_BufferResults( false ) -{ - this->m_Buffer.clear( ); -} - -// ------------------------------------------------------------------------- -template< class G > -cpExtensions::Algorithms::GradientFunctionBase< G >:: -~GradientFunctionBase( ) -{ - this->m_Buffer.clear( ); -} - -#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__HXX__ - -// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/GradientImageFunctionBase.h b/lib/cpExtensions/Algorithms/GradientImageFunctionBase.h new file mode 100644 index 0000000..d416e04 --- /dev/null +++ b/lib/cpExtensions/Algorithms/GradientImageFunctionBase.h @@ -0,0 +1,68 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__ +#define __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__ + +#include + +namespace cpExtensions +{ + namespace Algorithms + { + /** + * Base class to compute values based on image gradients (vector). + */ + template< class _TGradient > + class GradientImageFunctionBase + : public itk::ImageFunction< _TGradient, typename _TGradient::PixelType::ValueType, typename _TGradient::PixelType::ValueType > + { + public: + // Types from input arguments + typedef _TGradient TGradient; + typedef typename _TGradient::PixelType TVector; + typedef typename TVector::ValueType TScalar; + itkStaticConstMacro( Dimension, unsigned int, _TGradient::ImageDimension ); + + // Standard itk types + typedef GradientImageFunctionBase Self; + typedef itk::ImageFunction< _TGradient, TScalar, TScalar > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + // Types from base itk::ImageFunction + typedef typename Superclass::InputType TInput; + typedef typename Superclass::OutputType TOutput; + typedef typename Superclass::PointType TPoint; + typedef typename Superclass::ContinuousIndexType TContIndex; + typedef typename Superclass::IndexType TIndex; + + public: + itkTypeMacro( GradientImageFunctionBase, itkImageFunction ); + + public: + virtual void Prepare( ) const; + virtual TOutput Evaluate( const TPoint& p ) const ITK_OVERRIDE; + virtual TOutput EvaluateAtIndex( const TIndex& i ) const ITK_OVERRIDE; + virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const ITK_OVERRIDE; + + protected: + GradientImageFunctionBase( ); + virtual ~GradientImageFunctionBase( ); + + virtual TOutput _Evaluate( const TIndex& i ) const = 0; + + private: + // Purposely not implemented. + GradientImageFunctionBase( const Self& ); + void operator=( const Self& ); + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/GradientImageFunctionBase.hxx b/lib/cpExtensions/Algorithms/GradientImageFunctionBase.hxx new file mode 100644 index 0000000..f4cb928 --- /dev/null +++ b/lib/cpExtensions/Algorithms/GradientImageFunctionBase.hxx @@ -0,0 +1,59 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__HXX__ +#define __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__HXX__ + +// ------------------------------------------------------------------------- +template< class _TGradient > +void cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +Prepare( ) const +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +typename cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +TOutput cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +Evaluate( const TPoint& p ) const +{ + TIndex i; + this->GetInputImage( )->TransformPhysicalPointToIndex( p, i ); + return( this->EvaluateAtIndex( i ) ); +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +typename cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +TOutput cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +EvaluateAtIndex( const TIndex& i ) const +{ + return( this->_Evaluate( i ) ); +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +typename cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +TOutput cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +EvaluateAtContinuousIndex( const TContIndex& i ) const +{ + TPoint p; + this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p ); + return( this->Evaluate( p ) ); +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +GradientImageFunctionBase( ) + : Superclass( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::GradientImageFunctionBase< _TGradient >:: +~GradientImageFunctionBase( ) +{ +} + +#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__HXX__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/GulsunTekMedialness.h b/lib/cpExtensions/Algorithms/GulsunTekMedialness.h index 2fdf61c..49eb762 100644 --- a/lib/cpExtensions/Algorithms/GulsunTekMedialness.h +++ b/lib/cpExtensions/Algorithms/GulsunTekMedialness.h @@ -1,12 +1,7 @@ -// ------------------------------------------------------------------------- -// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) -// ------------------------------------------------------------------------- - #ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__ #define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__ -#include -#include +#include namespace cpExtensions { @@ -14,35 +9,31 @@ namespace cpExtensions { /** */ - template< class G > + template< class _TGradient > class GulsunTekMedialness - : public GradientFunctionBase< G > + : public GradientImageFunctionBase< _TGradient > { public: - // Standard itk types - typedef GulsunTekMedialness Self; - typedef GradientFunctionBase< G > Superclass; - typedef itk::SmartPointer< Self > Pointer; - typedef itk::SmartPointer< const Self > ConstPointer; - - // Types from superclass - typedef typename Superclass::TGradient TGradient; - typedef typename Superclass::TVector TVector; - typedef typename Superclass::TScalar TScalar; - typedef typename Superclass::TInput TInput; - typedef typename Superclass::TOutput TOutput; - typedef typename Superclass::TPoint TPoint; - typedef typename Superclass::TContIndex TContIndex; - typedef typename Superclass::TIndex TIndex; - typedef typename Superclass::TBuffer TBuffer; - - typedef typename TIndex::OffsetType TOffset; - typedef std::vector< double > TProfile; - typedef std::vector< TOffset > TOffsets; + typedef GulsunTekMedialness Self; + typedef GradientImageFunctionBase< _TGradient > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + itkStaticConstMacro( Dimension, unsigned int, Superclass::Dimension ); + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TScalar TScalar; + typedef typename Superclass::TIndex TIndex; + typedef typename Superclass::TVector TVector; + typedef typename Superclass::TPoint TPoint; + typedef typename TIndex::OffsetType TOffset; + + typedef std::vector< double > TProfile; + typedef std::vector< TOffset > TOffsets; public: itkNewMacro( Self ); - itkTypeMacro( GulsunTekMedialness, GradientFunctionBase ); + itkTypeMacro( GulsunTekMedialness, GradientImageFunctionBase ); itkGetConstMacro( MinRadius, double ); itkGetConstMacro( MaxRadius, double ); @@ -58,7 +49,7 @@ namespace cpExtensions GulsunTekMedialness( ); virtual ~GulsunTekMedialness( ); - virtual TOutput _Evaluate( const TIndex& i ) const; + virtual TOutput _Evaluate( const TIndex& i ) const ITK_OVERRIDE; private: // Purposely not implemented. @@ -77,7 +68,7 @@ namespace cpExtensions } // ecapseman #ifndef ITK_MANUAL_INSTANTIATION -#include +# include #endif // ITK_MANUAL_INSTANTIATION #endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__ diff --git a/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx b/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx index 01a278b..941a37c 100644 --- a/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx +++ b/lib/cpExtensions/Algorithms/GulsunTekMedialness.hxx @@ -1,24 +1,13 @@ -// ------------------------------------------------------------------------- -// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co) -// ------------------------------------------------------------------------- - #ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__ #define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__ #include #include - - -#include -#include -#include - - - +#include // ------------------------------------------------------------------------- -template< class G > -cpExtensions::Algorithms::GulsunTekMedialness< G >:: +template< class _TGradient > +cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >:: GulsunTekMedialness( ) : Superclass( ), m_MinRadius( double( 0 ) ), @@ -29,167 +18,95 @@ GulsunTekMedialness( ) } // ------------------------------------------------------------------------- -template< class G > -cpExtensions::Algorithms::GulsunTekMedialness< G >:: +template< class _TGradient > +cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >:: ~GulsunTekMedialness( ) { } // ------------------------------------------------------------------------- -template< class G > -typename cpExtensions::Algorithms::GulsunTekMedialness< G >:: -TOutput cpExtensions::Algorithms::GulsunTekMedialness< G >:: +template< class _TGradient > +typename cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >:: +TOutput cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >:: _Evaluate( const TIndex& i ) const { - const G* gradient = this->GetInputImage( ); - typename G::RegionType region = gradient->GetRequestedRegion( ); - typename G::PointType i_P; - gradient->TransformIndexToPhysicalPoint( i, i_P ); - - std::queue< TIndex > q; - std::set< TIndex, typename TIndex::LexicographicCompare > m; - std::map< TOffset, std::vector< TIndex >, typename TOffset::LexicographicCompare > profiles; - q.push( i ); - - while( !q.empty( ) ) + itk::Object::GlobalWarningDisplayOff( ); + + // Various values + const _TGradient* in = this->GetInputImage( ); + double pi2n = + double( 2 ) * double( vnl_math::pi ) / + double( this->m_ProfileSampling ); + double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 ); + double optR = double( 0 ); + TPoint pnt; + in->TransformIndexToPhysicalPoint( i, pnt ); + + // Main loop + for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ ) { - // Get next node - TIndex node = q.front( ); - q.pop( ); - if( m.find( node ) != m.end( ) ) - continue; - m.insert( node ); - - // Get neighborhood - for( unsigned int d = 0; d < Self::Dimension; d++ ) + for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ ) { - for( int off = -1; off <= 1; off += 2 ) + TProfile maxProfile( this->m_RadialSampling, double( 0 ) ); + for( unsigned int p = 0; p < this->m_ProfileSampling; p++ ) { - TIndex neighbor = node; - neighbor[ d ] += off; - if( region.IsInside( neighbor ) ) + double a = pi2n * double( p ); + + // Direction of this profile + TVector dir; + dir.Fill( TScalar( 0 ) ); + dir[ cx ] = TScalar( std::cos( a ) ); + dir[ cy ] = TScalar( std::sin( a ) ); + + double maxrise = double( 0 ); + double maxfall = double( -1 ); + TProfile profile; + for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) { - typename G::PointType neighbor_P; - gradient->TransformIndexToPhysicalPoint( neighbor, neighbor_P ); - if( double( i_P.EuclideanDistanceTo( neighbor_P ) ) > this->m_MaxRadius ) - continue; - - // Normalize offset in terms of a l1 (manhattan) distance - TOffset offset = neighbor - i; - // std::cout << "Offset: " << offset << std::endl; - - /* - int max_off = 0; - for( unsigned int o = 0; o < Self::Dimension; ++o ) - max_off += std::abs( offset[ o ] ); - if( max_off == 0 ) - continue; - bool normalized = true; - TOffset normalized_offset; - for( unsigned int o = 0; o < Self::Dimension && normalized; ++o ) - { - if( offset[ o ] % max_off == 0 ) - normalized_offset[ o ] = offset[ o ] / max_off; - else - normalized = false; - - } // rof - if( !normalized ) - normalized_offset = offset; - std::cout << "offset: " << normalized_offset << std::endl; - - // Update profiles - profiles[ normalized_offset ].push_back( neighbor ); - */ - - // Update queue - q.push( neighbor ); - - } // fi + double radius = double( r ) * rOff; + TIndex idx; + typename TPoint::VectorType aux; + aux.SetVnlVector( dir.GetVnlVector( ) ); + if( + in->TransformPhysicalPointToIndex( pnt + ( aux * radius ), idx ) + ) + { + TVector g = in->GetPixel( idx ); + double b = double( g.GetNorm( ) ); + if( double( g * dir ) < double( 0 ) ) + b *= double( -1 ); + maxrise = ( b > maxrise )? b: maxrise; + if( radius >= this->m_MinRadius ) + maxfall = ( b < maxfall )? b: maxfall; + profile.push_back( -b - maxrise ); + } + else + profile.push_back( double( 0 ) ); + + } // rof + + for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) + { + double E = profile[ r ] / -maxfall; + E = ( E < double( 0 ) )? double( 0 ): E; + E = ( E > double( 1 ) )? double( 1 ): E; + maxProfile[ r ] += E; + + } // rof + + } // rof + + for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) + { + double E = maxProfile[ r ] / double( this->m_RadialSampling ); + optR = ( E > optR )? E: optR; } // rof } // rof - } // elihw - // std::cout << "HOLA: " << profiles.size( ) << std::endl; - - - return( TOutput( 0 ) ); - /* TODO - const G* in = this->GetInputImage( ); - double pi2n = - double( 2 ) * double( vnl_math::pi ) / - double( this->m_ProfileSampling ); - double rOff = this->m_MaxRadius / double( this->m_RadialSampling - 1 ); - double optR = double( 0 ); - TPoint pnt; - in->TransformIndexToPhysicalPoint( i, pnt ); - - // Main loop - for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ ) - { - for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ ) - { - TProfile maxProfile( this->m_RadialSampling, double( 0 ) ); - for( unsigned int p = 0; p < this->m_ProfileSampling; p++ ) - { - double a = pi2n * double( p ); - - // Direction of this profile - TVector dir; - dir.Fill( TScalar( 0 ) ); - dir[ cx ] = TScalar( std::cos( a ) ); - dir[ cy ] = TScalar( std::sin( a ) ); - - double maxrise = double( 0 ); - double maxfall = double( -1 ); - TProfile profile; - for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) - { - double radius = double( r ) * rOff; - TIndex idx; - if( - in->TransformPhysicalPointToIndex( pnt + ( dir * radius ), idx ) - ) - { - TVector g = in->GetPixel( idx ); - double b = double( g.GetNorm( ) ); - if( double( g * dir ) < double( 0 ) ) - b *= double( -1 ); - maxrise = ( b > maxrise )? b: maxrise; - if( radius >= this->m_MinRadius ) - maxfall = ( b < maxfall )? b: maxfall; - profile.push_back( -b - maxrise ); - } - else - profile.push_back( double( 0 ) ); - - } // rof - - for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) - { - double E = profile[ r ] / -maxfall; - E = ( E < double( 0 ) )? double( 0 ): E; - E = ( E > double( 1 ) )? double( 1 ): E; - maxProfile[ r ] += E; - - } // rof - - } // rof - - for( unsigned int r = 0; r < this->m_RadialSampling; r++ ) - { - double E = maxProfile[ r ] / double( this->m_RadialSampling ); - optR = ( E > optR )? E: optR; - - } // rof - - } // rof - - } // rof - return( TScalar( optR ) ); - */ + } // rof + return( TScalar( optR ) ); } #endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__ diff --git a/lib/cpExtensions/Algorithms/MFluxMedialness.h b/lib/cpExtensions/Algorithms/MFluxMedialness.h new file mode 100644 index 0000000..e050392 --- /dev/null +++ b/lib/cpExtensions/Algorithms/MFluxMedialness.h @@ -0,0 +1,74 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__ +#define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__ + +#include + +namespace cpExtensions +{ + namespace Algorithms + { + /** + */ + template< class _TGradient > + class MFluxMedialness + : public GradientImageFunctionBase< _TGradient > + { + public: + typedef MFluxMedialness Self; + typedef GradientImageFunctionBase< _TGradient > Superclass; + typedef itk::SmartPointer< Self > Pointer; + typedef itk::SmartPointer< const Self > ConstPointer; + + itkStaticConstMacro( Dimension, unsigned int, Superclass::Dimension ); + + typedef typename Superclass::TOutput TOutput; + typedef typename Superclass::TScalar TScalar; + typedef typename Superclass::TIndex TIndex; + typedef typename Superclass::TVector TVector; + typedef typename Superclass::TPoint TPoint; + + typedef std::vector< double > TRCandidates; + + public: + itkNewMacro( Self ); + itkTypeMacro( MFluxMedialness, GradientImageFunctionBase ); + + itkGetConstMacro( RadiusStep, double ); + itkGetConstMacro( MinRadius, double ); + itkGetConstMacro( MaxRadius, double ); + itkGetConstMacro( RadialSampling, unsigned int ); + + itkSetMacro( RadiusStep, double ); + itkSetMacro( MinRadius, double ); + itkSetMacro( MaxRadius, double ); + itkSetMacro( RadialSampling, unsigned int ); + + protected: + MFluxMedialness( ); + virtual ~MFluxMedialness( ); + + virtual TOutput _Evaluate( const TIndex& i ) const ITK_OVERRIDE; + + private: + // Purposely not implemented. + MFluxMedialness( const Self& ); + void operator=( const Self& ); + + protected: + double m_MinRadius; + double m_MaxRadius; + unsigned int m_RadialSampling; + double m_RadiusStep; + }; + + } // ecapseman + +} // ecapseman + +#ifndef ITK_MANUAL_INSTANTIATION +# include +#endif // ITK_MANUAL_INSTANTIATION + +#endif // __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__ + +// eof - $RCSfile$ diff --git a/lib/cpExtensions/Algorithms/MFluxMedialness.hxx b/lib/cpExtensions/Algorithms/MFluxMedialness.hxx new file mode 100644 index 0000000..20e7f15 --- /dev/null +++ b/lib/cpExtensions/Algorithms/MFluxMedialness.hxx @@ -0,0 +1,159 @@ +#ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__ +#define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__ + +#include +#include +#include + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::MFluxMedialness< _TGradient >:: +MFluxMedialness( ) + : Superclass( ), + m_MinRadius( double( 0 ) ), + m_MaxRadius( double( 1 ) ), + m_RadialSampling( 4 ), + m_RadiusStep( double( 1 ) ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +cpExtensions::Algorithms::MFluxMedialness< _TGradient >:: +~MFluxMedialness( ) +{ +} + +// ------------------------------------------------------------------------- +template< class _TGradient > +typename cpExtensions::Algorithms::MFluxMedialness< _TGradient >:: +TOutput cpExtensions::Algorithms::MFluxMedialness< _TGradient >:: +_Evaluate( const TIndex& i ) const +{ + itk::Object::GlobalWarningDisplayOff( ); + + double pi2n = double( 2 ) * double( vnl_math::pi ); + pi2n /= int( this->m_RadialSampling ); + const _TGradient* img = this->GetInputImage( ); + //const itk::Image::SpacingType& input_spacing = img->GetSpacing( ); + + double Flux1 = 0; + double Flux2 = 0; + double MFlux = 0; + + TRCandidates FluxFinal; + TRCandidates radiusGenerated; + double dR = double( 0 ); + double optR = double( 0 ); + TPoint center; + img->TransformIndexToPhysicalPoint( i, center ); + double radius; + + for( unsigned int cx = 0; cx < Self::Dimension - 1; cx++ ) + { + for( unsigned int cy = cx + 1; cy < Self::Dimension; cy++ ) + { + dR = double( 0 ); + FluxFinal.clear(); + radiusGenerated.clear(); + radius = this->m_MinRadius; + while( radius <= this->m_MaxRadius ) + { + MFlux = 0; + for( unsigned int I_radial = 0; I_radial < this->m_RadialSampling / 2; I_radial++ ) + { + Flux1 = 0; + Flux2 = 0; + + // Direction of first profile + typename TPoint::VectorType dir1; + dir1.Fill( double( 0 ) ); + dir1[ cx ] = std::cos( pi2n * double( I_radial ) ); + dir1[ cy ] = std::sin( pi2n * double( I_radial ) ); + //dir1 *= (radius); + + TIndex rIdx; + + if ( img->TransformPhysicalPointToIndex( center + (dir1*radius), rIdx ) ) + { + TVector grad_rIdx = img->GetPixel( rIdx ); + TVector u_i1; + u_i1.SetVnlVector( ( center - ( center + dir1 ) ).GetVnlVector( ) ); + u_i1.Normalize( ); + // dot product + Flux1 = grad_rIdx * u_i1; + } + else + { + //if (Self::Dimension==3) + //{ + //std::cout<<"Point Edge x:"<TransformPhysicalPointToIndex( center + (dir2*radius), rIdx2 ) ) + { + TVector grad_rIdx2 = img->GetPixel( rIdx2 ); + TVector u_i2; + u_i2.SetVnlVector( ( center - ( center + dir2 ) ).GetVnlVector( ) ); + u_i2.Normalize( ); + + Flux2 = grad_rIdx2 * u_i2; + } + else + { + //if (Self::Dimension==3) + //{ + //std::cout<<"Point Edge x:"<m_RadialSampling; + FluxFinal.push_back(MFlux); + radiusGenerated.push_back(radius); + + radius += this->m_RadiusStep; + + } //elihw + + dR= *( std::max_element( FluxFinal.begin(), FluxFinal.end() ) ); + optR= (dR>optR)? dR:optR; + + } // rof + + } // rof + return( TScalar(optR) ); +} + +#endif // __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__ + +// eof - $RCSfile$ -- 2.47.1