#include <itkImageRegionConstIteratorWithIndex.h>
namespace cpExtensions
+{
+ namespace Algorithms
{
- namespace Algorithms
+ /**
+ */
+ template< class I, class S = float >
+ class InertiaMedialness
+ : public itk::ImageFunction< I, S, S >
{
- /**
- */
- template< class I, class S = float >
- class InertiaMedialness
- : public itk::ImageFunction< I, S, S >
- {
- public:
- // Standard itk types
- typedef InertiaMedialness Self;
- typedef itk::ImageFunction< I, S, S > 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;
- typedef typename TIndex::OffsetType TOffset;
-
- // Sparse buffer
- typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer;
-
- typedef InertiaTensorFunction< S, I::ImageDimension > TInertia;
-
- public:
- itkNewMacro( Self );
- itkTypeMacro( InertiaMedialness, itkImageFunction );
-
- itkBooleanMacro( BufferResults );
- itkGetConstMacro( BufferResults, bool );
- itkGetConstMacro( MaxRadius, double );
-
- itkSetMacro( BufferResults, bool );
- itkSetMacro( MaxRadius, double );
-
- public:
- virtual void ResetBuffer( )
+ public:
+ // Standard itk types
+ typedef InertiaMedialness Self;
+ typedef itk::ImageFunction< I, S, S > 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;
+ typedef typename TIndex::OffsetType TOffset;
+
+ // Sparse buffer
+ typedef std::map< TIndex, TOutput, typename TIndex::LexicographicCompare > TBuffer;
+
+ typedef InertiaTensorFunction< S, I::ImageDimension > TInertia;
+
+ public:
+ itkNewMacro( Self );
+ itkTypeMacro( InertiaMedialness, itkImageFunction );
+
+ itkBooleanMacro( BufferResults );
+ itkGetConstMacro( BufferResults, bool );
+ itkGetConstMacro( MaxRadius, double );
+
+ itkSetMacro( BufferResults, bool );
+ itkSetMacro( MaxRadius, double );
+
+ public:
+ virtual void ResetBuffer( )
+ {
+ this->m_Buffer.clear( );
+ }
+
+ virtual TOutput Evaluate( const TPoint& p ) const
+ {
+ TIndex i;
+ this->GetInputImage( )->TransformPhysicalPointToIndex( p, i );
+ return( this->EvaluateAtIndex( i ) );
+ }
+
+ virtual TOutput EvaluateAtIndex( const TIndex& i ) const
+ {
+ TOutput res = TOutput( 0 );
+ bool computed = false;
+ if( this->m_BufferResults )
{
- this->m_Buffer.clear( );
- }
-
- virtual TOutput Evaluate( const TPoint& p ) const
+ 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 );
+ }
+
+ virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const
+ {
+ TPoint p;
+ this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p );
+ return( this->Evaluate( p ) );
+ }
+
+ protected:
+ InertiaMedialness( )
+ : Superclass( ),
+ m_BufferResults( false ),
+ m_MaxRadius( double( 1 ) )
+ {
+ this->m_Buffer.clear( );
+ }
+
+ virtual ~InertiaMedialness( )
+ {
+ this->m_Buffer.clear( );
+ }
+
+ virtual TOutput _Evaluate( const TIndex& idx ) const
+ {
+ const I* image = this->GetInputImage( );
+
+ typename I::PointType p_i;
+ image->TransformIndexToPhysicalPoint( idx, p_i );
+
+ typename I::PointType max_p, min_p;
+ for( unsigned int d = 0; d < I::ImageDimension; ++d )
{
- TIndex i;
- this->GetInputImage( )->TransformPhysicalPointToIndex( p, i );
- return( this->EvaluateAtIndex( i ) );
- }
-
- virtual TOutput EvaluateAtIndex( const TIndex& i ) const
+ max_p[ d ] = p_i[ d ] + this->m_MaxRadius;
+ min_p[ d ] = p_i[ d ] - this->m_MaxRadius;
+
+ } // rof
+ TIndex max_i, min_i;
+ image->TransformPhysicalPointToIndex( max_p, max_i );
+ image->TransformPhysicalPointToIndex( min_p, min_i );
+
+ typename I::RegionType in_region = image->GetRequestedRegion( );
+ TIndex in_index = in_region.GetIndex( );
+ TIndex in_last = in_index + in_region.GetSize( );
+ typename I::SizeType size;
+ for( unsigned int d = 0; d < I::ImageDimension; ++d )
{
- 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;
+ if( min_i[ d ] < in_index[ d ] ) min_i[ d ] = in_index[ d ];
+ if( max_i[ d ] < in_index[ d ] ) max_i[ d ] = in_index[ d ];
+ if( min_i[ d ] >= in_last[ d ] ) min_i[ d ] = in_last[ d ];
+ if( max_i[ d ] >= in_last[ d ] ) max_i[ d ] = in_last[ d ];
- } // fi
+ size[ d ] = max_i[ d ] - min_i[ d ];
- if( !computed )
- res = this->_Evaluate( i );
+ } // rof
- if( this->m_BufferResults )
- this->m_Buffer[ i ] = res;
- return( res );
- }
+ typename I::RegionType region;
+ region.SetIndex( min_i );
+ region.SetSize( size );
- virtual TOutput EvaluateAtContinuousIndex( const TContIndex& i ) const
+ std::vector< typename TInertia::Pointer > inertias;
+ itk::ImageRegionConstIteratorWithIndex< I > it( image, region );
+ for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
{
- TPoint p;
- this->GetInputImage( )->TransformContinuousIndexToPhysicalPoint( i, p );
- return( this->Evaluate( p ) );
- }
+ TOffset off = it.GetIndex( ) - idx;
+ unsigned long l1dist = std::abs( off[ 0 ] );
+ for( unsigned int d = 1; d < I::ImageDimension; ++d )
+ l1dist = ( std::abs( off[ d ] ) > l1dist )? std::abs( off[ d ] ): l1dist;
- protected:
- InertiaMedialness( )
- : Superclass( ),
- m_BufferResults( false ),
- m_MaxRadius( double( 1 ) )
- {
- this->m_Buffer.clear( );
- }
+ typename TInertia::TPoint i_pnt;
+ image->TransformIndexToPhysicalPoint( it.GetIndex( ), i_pnt );
- virtual ~InertiaMedialness( )
- {
- this->m_Buffer.clear( );
- }
-
- virtual TOutput _Evaluate( const TIndex& idx ) const
- {
- const I* image = this->GetInputImage( );
-
- typename I::PointType p_i;
- image->TransformIndexToPhysicalPoint( idx, p_i );
-
- typename I::PointType max_p, min_p;
- for( unsigned int d = 0; d < I::ImageDimension; ++d )
+ for( unsigned long l = 0; l < l1dist; ++l )
{
- max_p[ d ] = p_i[ d ] + this->m_MaxRadius;
- min_p[ d ] = p_i[ d ] - this->m_MaxRadius;
+ if( inertias.size( ) <= l )
+ inertias.push_back( TInertia::New( ) );
+ inertias[ l ]->AddMass( i_pnt.GetVectorFromOrigin( ), S( it.Get( ) ) );
+
+ /* TODO
+ typename TInertias::iterator inIt = inertias.find( l );
+ if( inIt == inertias.end( ) )
+ inIt = inertias.insert( std::pair< unsigned long, typename TInertia::Pointer >( l, TInertia::New( ) ) ).first;
+ */
} // rof
- TIndex max_i, min_i;
- image->TransformPhysicalPointToIndex( max_p, max_i );
- image->TransformPhysicalPointToIndex( min_p, min_i );
-
- typename I::RegionType in_region = image->GetRequestedRegion( );
- TIndex in_index = in_region.GetIndex( );
- TIndex in_last = in_index + in_region.GetSize( );
- typename I::SizeType size;
- for( unsigned int d = 0; d < I::ImageDimension; ++d )
- {
- if( min_i[ d ] < in_index[ d ] ) min_i[ d ] = in_index[ d ];
- if( max_i[ d ] < in_index[ d ] ) max_i[ d ] = in_index[ d ];
- if( min_i[ d ] >= in_last[ d ] ) min_i[ d ] = in_last[ d ];
- if( max_i[ d ] >= in_last[ d ] ) max_i[ d ] = in_last[ d ];
- size[ d ] = max_i[ d ] - min_i[ d ];
+ } // rof
- } // rof
-
- typename I::RegionType region;
- region.SetIndex( min_i );
- region.SetSize( size );
-
- std::vector< typename TInertia::Pointer > inertias;
- itk::ImageRegionConstIteratorWithIndex< I > it( image, region );
- for( it.GoToBegin( ); !it.IsAtEnd( ); ++it )
+ if( inertias.size( ) > 0 )
+ {
+ S res = S( 0 );
+ for( unsigned int l = 0; l < inertias.size( ); ++l )
{
- TOffset off = it.GetIndex( ) - idx;
- unsigned long l1dist = std::abs( off[ 0 ] );
- for( unsigned int d = 1; d < I::ImageDimension; ++d )
- l1dist = ( std::abs( off[ d ] ) > l1dist )? std::abs( off[ d ] ): l1dist;
-
- typename TInertia::TPoint i_pnt;
- image->TransformIndexToPhysicalPoint( it.GetIndex( ), i_pnt );
-
- for( unsigned long l = 0; l < l1dist; ++l )
- {
- if( inertias.size( ) <= l )
- inertias.push_back( TInertia::New( ) );
- inertias[ l ]->AddMass( i_pnt.GetVectorFromOrigin( ), S( it.Get( ) ) );
-
- /* TODO
- typename TInertias::iterator inIt = inertias.find( l );
- if( inIt == inertias.end( ) )
- inIt = inertias.insert( std::pair< unsigned long, typename TInertia::Pointer >( l, TInertia::New( ) ) ).first;
- */
-
- } // rof
+ typename TInertia::TVector pv, r;
+ typename TInertia::TMatrix pm;
+ inertias[ l ]->GetEigenAnalysis( pm, pv, r );
+ S v = pv.GetNorm( );
+ if( l == 0 || v > res )
+ res = v;
} // rof
-
- if( inertias.size( ) > 0 )
- {
- S res = S( 0 );
- for( unsigned int l = 0; l < inertias.size( ); ++l )
- {
- typename TInertia::TVector pv, r;
- typename TInertia::TMatrix pm;
- inertias[ l ]->GetEigenAnalysis( pm, pv, r );
- S v = pv.GetNorm( );
- if( l == 0 || v > res )
- res = v;
-
- } // rof
- return( res );
- }
- else
- return( TOutput( 0 ) );
+ return( res );
}
+ else
+ return( TOutput( 0 ) );
+ }
- private:
- // Purposely not implemented.
- InertiaMedialness( const Self& );
- void operator=( const Self& );
+ private:
+ // Purposely not implemented.
+ InertiaMedialness( const Self& );
+ void operator=( const Self& );
- protected:
- mutable TBuffer m_Buffer;
- bool m_BufferResults;
+ protected:
+ mutable TBuffer m_Buffer;
+ bool m_BufferResults;
- double m_MaxRadius;
- };
+ double m_MaxRadius;
+ };
- } // ecapseman
+ } // ecapseman
} // ecapseman
+#ifndef ITK_MANUAL_INSTANTIATION
// TODO: #include <cpExtensions/Algorithms/InertiaMedialness.hxx>
+#endif // ITK_MANUAL_INSTANTIATION
#endif // __CPEXTENSIONS__ALGORITHMS__INERTIAMEDIALNESS__H__
git://git.creatis.insa-lyon.fr / cpPlugins.git / blobdiff