--- /dev/null
+#ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__
+#define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__
+
+#include <cpExtensions/Algorithms/GradientImageFunctionBase.h>
+
+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 <cpExtensions/Algorithms/FluxMedialness.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__H__
+
+// eof - $RCSfile$
--- /dev/null
+#ifndef __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__
+#define __CPEXTENSIONS__ALGORITHMS__FLUXMEDIALNESS__HXX__
+
+#include <cmath>
+#include <vnl/vnl_math.h>
+#include <itkLineConstIterator.h>
+
+// -------------------------------------------------------------------------
+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:"<<radius<<std::endl;
+ //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
+ //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
+ Flux /= this->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$
+++ /dev/null
-// -------------------------------------------------------------------------
-// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
-// -------------------------------------------------------------------------
-
-#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__
-#define __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__
-
-#include <map>
-#include <itkImageFunction.h>
-
-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 <cpExtensions/Algorithms/GradientFunctionBase.hxx>
-#endif // ITK_MANUAL_INSTANTIATION
-
-#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTFUNCTIONBASE__H__
-
-// eof - $RCSfile$
+++ /dev/null
-// -------------------------------------------------------------------------
-// @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$
--- /dev/null
+#ifndef __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__
+#define __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__
+
+#include <itkImageFunction.h>
+
+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 <cpExtensions/Algorithms/GradientImageFunctionBase.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __CPEXTENSIONS__ALGORITHMS__GRADIENTIMAGEFUNCTIONBASE__H__
+
+// eof - $RCSfile$
--- /dev/null
+#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$
-// -------------------------------------------------------------------------
-// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
-// -------------------------------------------------------------------------
-
#ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__
#define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__
-#include <vector>
-#include <cpExtensions/Algorithms/GradientFunctionBase.h>
+#include <cpExtensions/Algorithms/GradientImageFunctionBase.h>
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 );
GulsunTekMedialness( );
virtual ~GulsunTekMedialness( );
- virtual TOutput _Evaluate( const TIndex& i ) const;
+ virtual TOutput _Evaluate( const TIndex& i ) const ITK_OVERRIDE;
private:
// Purposely not implemented.
} // ecapseman
#ifndef ITK_MANUAL_INSTANTIATION
-#include <cpExtensions/Algorithms/GulsunTekMedialness.hxx>
+# include <cpExtensions/Algorithms/GulsunTekMedialness.hxx>
#endif // ITK_MANUAL_INSTANTIATION
#endif // __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__H__
-// -------------------------------------------------------------------------
-// @author Leonardo Florez-Valencia (florez-l@javeriana.edu.co)
-// -------------------------------------------------------------------------
-
#ifndef __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
#define __CPEXTENSIONS__ALGORITHMS__GULSUNTEKMEDIALNESS__HXX__
#include <cmath>
#include <vnl/vnl_math.h>
-
-
-#include <queue>
-#include <map>
-#include <set>
-
-
-
+#include <itkLineConstIterator.h>
// -------------------------------------------------------------------------
-template< class G >
-cpExtensions::Algorithms::GulsunTekMedialness< G >::
+template< class _TGradient >
+cpExtensions::Algorithms::GulsunTekMedialness< _TGradient >::
GulsunTekMedialness( )
: Superclass( ),
m_MinRadius( double( 0 ) ),
}
// -------------------------------------------------------------------------
-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__
--- /dev/null
+#ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__
+#define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__
+
+#include <cpExtensions/Algorithms/GradientImageFunctionBase.h>
+
+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 <cpExtensions/Algorithms/MFluxMedialness.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
+
+#endif // __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__H__
+
+// eof - $RCSfile$
--- /dev/null
+#ifndef __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
+#define __CPEXTENSIONS__ALGORITHMS__MFLUXMEDIALNESS__HXX__
+
+#include <cmath>
+#include <vnl/vnl_math.h>
+#include <itkLineConstIterator.h>
+
+// -------------------------------------------------------------------------
+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:"<<center[0]+dir1[0] ;
+ //std::cout<<" y:"<<center[1]+dir1[1]<<" z:"<<center[2]+dir1[2]<<std::endl;
+ //}
+ //else if (Self::Dimension==2)
+ //{
+ //std::cout<<"Point Edge x:"<<center[0]+dir1[0] ;
+ //std::cout<<" y:"<<center[1]+dir1[1]<<std::endl;
+ //}
+ }
+
+ // Direction of second profile
+ // pi2n*Iradial + 180°
+ typename TPoint::VectorType dir2;
+ dir2.Fill( double( 0 ) );
+ dir2[ cx ] = std::cos( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
+ dir2[ cy ] = std::sin( (pi2n) * double( I_radial ) + double( vnl_math::pi ));
+
+ TIndex rIdx2;
+
+ if ( img->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:"<<center[0]+dir2[0] ;
+ //std::cout<<" y:"<<center[1]+dir2[1]<<" z:"<<center[2]+dir2[2]<<std::endl;
+ //}
+ //else if (Self::Dimension==2)
+ //{
+ //std::cout<<"Point Edge x:"<<center[0]+dir2[0] ;
+ //std::cout<<" y:"<<center[1]+dir2[1]<<std::endl;
+ //}
+ }
+
+ MFlux += std::min( Flux1, Flux2 );
+ } // rof
+
+ //std::cout<<Self::Dimension<<" radius:"<<radius<<std::endl;
+ //std::cout<<"Center:"<<center[0]<<" "<<center[1]<<std::endl;
+ //std::cout<<"edge:"<<center[0]+radius*std::cos( pi2n * double( 0 ) )<<std::endl;
+
+ MFlux *= 2;
+ MFlux /= this->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$
git://git.creatis.insa-lyon.fr / cpPlugins.git / commitdiff